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CALCULATION  OF  THE  LIFT  DISTRIBUTION  AND  AERODYNAMIC  DERIVATIVES 
OF  QUASI-STATIC  ELASTIC  AIRCRAFT 


Liu  Qiangang,  Wu  Changlin»  Jian  Zheng 
Northwestern  Polytechnical  University 

Abstract : 

A  numerical  metbodHb  1>i€3eted  for  predicting  the  aerodynamic  characte¬ 
ristics  of  elastic  aircraft  under  quasi-static  approximation.  This  method  can  be 
used  to  eraluate  the  lift  distribution  and  11  main  longitudinal  aerodynamic 
derivatives  of  elastic  aircraft  at  subsonic  speeds.  The  aerodynamic  calculations 
are  based  on  the  Green's  function  method.  The  structure  deformation  U  evalua¬ 
ted  bp  the  free  structure  influence  coefficient  method.  The  cor.uination  of 
tbMe  methods  can  provide  an  efficient,  general  and  flexible  aerodynamic  tool 
for  design  of  elastic  aircraft. 

Several  numerical  examples  are  given  and  some  dynamical  problems  of  elastic 
aircraft  are  also  discussed  in.  this  paper.  The  derivatives  evaluated  in  this  paper 
can  be  directly  adopted  ifl'atialysis  of  stability  and  control  of  elastic  aircraft 

I*  FOREWORD 

Numerical  value  calculation  method  is  one  important  means  of  the 
present  research  on  elastic  aircraft  pneumatic  characteristics.  Be¬ 
cause  this  calculation  involves  problems  in  the  areas  of  pneumatics, 
structure,  and  flight  mechanics,  the  amount  of  work  is  comparatively 
large,  and  we  must  choose  the  appropriate  method. 

In  the  area  of  pneumatic  force  calculations,  this  text  uses  the 
speed  force  integral  equation  method  that  Greene’s  theorem  gives  guid¬ 
ance  to.  Its  major  merits  are  that  it  can  be  applied  to  pneumatic  force 
computations  of  various  complex  aircraft  contours,  and  excess  machine 
time  is  comparatively  less,  and  the  results  are  relatively  accurate. 

Thus  it  is  comparatively  appropriate  for  pneumatic  force  computations 
of  elastic  flight  instruments.  In  the  area  of  structural  deformation 
computations,  we  used  the  free  structure  influence  coefficient  method. 
This  method  not  only  correctly  explains  the  computation  problems  of 
aircraft  structure  deformations,  but  also  because  of  the  complete  a- 
likeness  of  the  elastic  aircraft  pneumatic  derivative  that  is  gained 


1 


f 


ca 


and  other  applicable  dynamic  equations  on  the  form  it  analyzes  and 
researches  the  dynamic  state  characteristics  of  elastic  aircraft  to 
bring  many  conveniences. 

In  this  text  we  will  carry  out  research  and  discussion  on  the 
strong  question  of  when  to  use  the  influence  coefficient  method  to 
calculate  aircraft  structure  deformation,  the  problem  of  the  co¬ 
ordinate  axis  system,  and  existing  problems  in  some  related  documents 
of  foreign  research  of  elastic  aircraft  dynamic  state  characteristics. 

II.  TURBULENCE  SPEED  FORCE  EQUATIONS  AND  OTHER  NUMERIC  VALUE 
EVALUATIONS 


FIG.  1:  COORDINATE  SYSTEM 


In  pneumatic  calculations,  assume  the  coordinate  system  in  Fig.  1 
(in  which  ox  axis  and  undisturbed  airflow  velocity  directions  are 
the  same). 


By  the  literature  £l,2,  and  3}  when  airflow  flows  through  the 
aircraft,  on  the  basis  of  Green's  theorem,  on  the  aircraft  surface 
disturbance  speed  force  $  satisfies  the  following  integral  equa- 

k.  *>-- J J sip/ J • 


s,  s. 


(1) 


Sm 


r* 


Where 


X-s/0i,  K-*/*,  Z-*/5,  0-  ✓ 1  -ML 


<*>  is  the  disturbance  speed  for?®  is  the  disturbance  speed  force 
method  toward  the  derivative;  c  is  the  average  pneumatic  chord  length; 
S0  expresses  the  aircraft  surface;  is  the  turbulence  area  that  is 
eased  out  by  the  wing  or  tail's  rear  edge;f  ,n,  and£are  the  integral 
variable  elements  of  the  directions  X,  Y,  and  Z  edges. 

To  design  the  aircraft  surface  equations  we  can  use: 


S(X,  Y,  Z)-  0  ( 

to  express,  thus  on  S  we  have  the  following  boundary  conditions: 


_! _  *s  i 

W3T  »x  ~fr 


(3) 


When  using  the  numeric  value  solution  to  explain  equation  (l), 
with  the  aircraft  surface  divides  into  N  small  4  side  form  elements, 
consider  the  speed  force  4>  value  on  every  element  as  equal  to  constant, 
and  also  equal  to  the  above-mentioned value  of  the  element's  center. 
On  every  element,  the  value  of  is  also  equal  to  constant,  and  also 
equal  to  the  value  of  the  above-mentioned  center  of  the  element 
The  vortex  area  divides  into  some  parallel  vortex  strips  with  the 
OX  axis.  Because  the  vortex  surface  cannot  bear  pressure,  each  vor¬ 
tex  strip  has: 


constant 


(4) 


Wh  ©  x*  ©  A 4* 

expresses  the  difference  of  the  vortex  strips  top  and  bottom 
surface  speed  force;  A$»«  expresses  the  difference  of  wing  surface 
rear  edge  top  and  bottom  surface  speed  force.  At  the  same  time,  on  the 
basis  of  storage  tower  conditions,  we  can  similarly  consider  that  A«&w 
is  equal  to  the  difference  of  the  element  centroid  top  and  bottom 
surface  speed  force  that  are  linked  together  with  the  above-mentioneda 
rear  edge,  to  ensure  that  the  rear  edge  pressure  difference  is  zero. 

Through  the  preceding  treatment  eaquation  (1)  can  turn  into  an 
N  exponential  linear  algebra  equation  that  possesses  N  unknown  mass$»: 


(5) 


Where 

&</  =  0. 


in,  is  an  enumerable  function.  When  i=j.  = 
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When  i4j, 

(6) 

(7) 

(8) 

(9) 


Where  R  is  the  distance  of  the  required  integral  point  by  the 

centroid  of  element  i;  Sg .  expresses  the  surface  of  element  j; 
expresses  the  vortex  strip  surface  of  the  rear  of  element  j. 


After  given  matter  surface  geometric  parameters,  we  can  obtain 


the  value  of  coefficients  b 


•  • ,  v_>  •  •  ■  W .  .  , 

ij  -iJ  iJ 


and  by  the  literature 
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j.  By  substituting  them  -in  equation  ^5)  and  evaluating  this  .equa¬ 
tion,  we  thus  obtain  the  value  $f(j=1,  2,  **’W)  of  each  element's 
speed  force.  Later  according  to  the  following  formula  we  can  obtain 
the  pressure  coefficient  of  each  element: 


*X>\ 


do) 


Where 


is  obtained  by  use  of  the  limited  deviation  method. 


After  obtaining  the  pressure  coefficient,  the  lift  and  pitching 
movement  coefficient  are  obtained  by  using  the  numeric  value  integral 
method.  The  deduction  resistance  coefficient  uses  the  attainment  of 
the  "joint  flow  *ield"  method.  These  coefficient  derivatives  can  be 
obtained  by  using  the  method  of  calculating  the  increments  of  these 
coefficients  and  the  ratio  of  corresponding  conditional  derivative 
increments . 


As  to  the  0  derivative,  it  also  can  be  obtained  under  the  con¬ 
dition  of  determined  common  assumptions.  To  design  the  center  of  mass 
of  the  aircraft  coil,  we  use  the  6  angle  velocity  to  make  the  pitching 
rotations,  thus  the  attached  angle  of  attack  of  the  centroid  of  ele¬ 
ment  j  that  has  emerged  is: 
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(11) 


This  condition  is  equivalent  to  the  aircraft  not  rotating  and  the 
dispersed  variable  curves  of  its  geometric  contour  according  to  . 
To  the  corresponding  curved  aircraft,  the  boundary  conditions  cahnge 
to : 

*H  -  -  (W.<co6Aa,+ r.jsta A«,)/0  (12) 

Where  Nxj  and  Txj  separately  express  the  normal  line  of  element  j 
and  the  cosine  of  the  tangent’s  same  airflow  direction  and  the  in¬ 
cluded  angle  of  the  OX  axis. 


With  formula  (12)  substiting  equation  (5),  later  according  to 
the  preceding  same  kind  of  measure  to  seek  explanation  we  can  obtain 
the  pneumatic  coefficient  when  the  pitching  angle  velocity  increments 
are  0,  and  with  the  mutual  elimination  of  this  pneumatic  coefficient 
increment  separation  and  the  dimensionless  angle  velocity  increments,* 
we  then  can  obtain  the  corresponding  pneumatic  derivatives. 


III.  SOME  RELATED  PROBLEMS  OF  COMPUTING  FREE  IN  FLIGHT  AIRCRAFT 
STRUCTURE  ELASTIC  DEFORMATION 


1 .  Quasi-Static  Assumptions  and  Other  Conditions 

When  researching  elastic  aircraft  dynamic  state  characteristics, 
in  order  to  simplify  computations  we  can  adopt  quasi-static  assumptions. 
So  called  quasi-static  assumptions  refer  to  when  researching  elastic 
aircraft  dynamic  state  characteristics , we  omit  a  degree  of  free  air¬ 
craft  elastic  movement,  and  consider  each  instantaneous  elastic  de¬ 
formation  of  aircraft  only  related  to  the  preceding  instantaneous 
effect  on  the  exterior  of  the  aircraft.  Also,  consider  that  there  is 
no  mutual  position  difference  between  force  and  deformation.  Except 
for  wanting  to  consider  the  special  characteristics  of  free  flight, 
the  relationship  between  them  can  be  completely  in  accordance  with  the 
treatment  of  the  relation  between  general  static  elastic  mechanics  cen¬ 
ter  force  and  deformation. 
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Adopting  quasi-static  assumptions  can  greatly  simplify  compu- 
taion.  But  only  if  the  free  oscillating  frequency  of  the  aircraft 
structure  is  comparatively  higher  than  the  rigid  body  motion  freq¬ 
uency  of  the  aircraft  (for  example  the  frequency  of  the  aircraft 
vertical  short  cycle  motion),  then  can  they  be  applied.  Because  at 
this  moment  the  coupling  between  the  elastic  motion  of  the  aircraft 
structure  and  the  rigid  body  motion  of  the  aircraft  is  very  weak.  When 
the  preceding  two  kinds  of  motion  frequency  are  comparatively  close, 
if  we  adopt  quasi-static  assumptions  then  we  possibly  will  have  fairly 
large  error.  At  this  time  to  research  the  dynamic  state  characterist¬ 
ics  of  aircraft  we  should  compute  the  degree  of  free  aircraft  elastic 
motion. 

2.  The  Average  Axis  System 

When  computing  free  flight  aircraft  structure  elastic  deformation 
this  text  uses  the  average  axis  system,  now  introduced  below.  When 
the  balanced  state  particle  system  of  the  free  motion  center  receives 
the  attached  exterior  effect,  its  disturbance  motion  can  be  broken 
down  into  three  sections:  a)  Each  section  of  the  particle  system  with 
its  center  of  mass  together  make  attached  translational  motion,  b) 

Each  section  of  the  particle  system  wind  its  center  of  mass  to  make 
rotating  motion,  c)  the  deformation  motion  of  the  particle  system. 

The  first  two  kinds  of  motion  are  jointly  called  rigid  body  motion. 

In  these  two  kinds  of  motion,  the  relative  position  between  each  par¬ 
ticle  of  the  particle  system  doesn’t  change.  The  third  kind  of  motion 
is  elastic  motion.  In  this  kind  of  motion,  the  relative  postion  of 
each  particle  changes.  Designing  5>i  to  separately  express  the 

displacement  of  certain  paticles  in  the  particle  system  in  translation 
al,  rotational,  and  deformation  motion.  Thus  the  total  displacement 
of  certain  particles  in  the  particle  system  is  equal  to: 


5<“5«+8*f+5«# 


(13) 


Although  the  relative  position  of  each  particle  within  the  par- 
tical  system  in  deformation  motion  do  change,  this  kind  of  change 
should  not  lead  to  the  displacement  of  the  position  of  the  center  of 
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mass  of  the  particle  system,  and  also  should  not  lead  to  the  whole 
particle  system  winding  its  center  of  mass  to  produce  rotational 
motion.  Because  these  two  kinds  of  motion  already  break  down  among 
the  former  two  kinds  of  rigid  body  motion.  Because  rigid  body  motion 
model  state  and  elastic  motion  model  state  are  mutually  perpendicular, 
the  elastic  displacement  of  each  particle  in  the  particle  system  5*« 
should  satisfy  the  following  two  conditions: 

N 

2  m, 8*i  “  0 

i  -  1 

N 

«-i 

Where  nu  is  the  mass  of  the  i  particle;  is  the  direct  loss  oi  .a 

i  particle  to  the  center  mass  of  the  particle  system;  N  is  the  s 
total  of  the  particles  in  the  particle  system. 

Now  we  assume  on  the  particle  system  selecting  this  kind  of 
coordinate  system  oxyz  with  one  right  angle,  when  the  particle  system 
does  not  receive  disturbance,  its  primary  position  is  on  the  center 
mass  of  the  particle  system,  and  the  direction  of  its  three  coordinate 
axis  is  determined  by  the  conditions  of  (15).  In  the  process  of 
disturbance  motion,  it  together  with  the  particle  system  makes  rigid 
body  motion.  This  is  also  to  say,  its  primary  point  throughout  mutually 
coincides  with  the  particle  system's  instantaneous  center  of  mass. 
Moreover,  it  rotates  together  with  the  aprticle  system.  Clearly,  this 
coordinate  system  is  a  non-inertia  coordinate  system.  At  the  same 
time  we  say  to  this  coordinate  system,  the  rigid  body  displacement  of 
each  particle  in  the  particle  system  of  course  can  only  be  zero. 
Therefore  the  displacement  of  each  particle  in  the  particle  system 
facing  this  coordinate  system  then  is  the  elastic  displacement  fa. 
that  we  require. 

This  kind  of  coordinate  system  is  called  the  average  body  axis 
system.  Its  position  can  be  determined  according  to  conditions  (14) 
and  (15). 

It  is  worthy  to  point  out  that  when  computing  structure  displace- 


(14) 
(1  ' 
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mer.t  the  structure  influence  coefficient  that  is  used  refers  to  its 
axis  system  as  the  structure  axis  system  (  its  coordinate  primary 
point  and  certain  structure  points  of  the  aircraft  are  consolidated). 
It  is  different  from  the  average  body  axis  system. 

IV.  COMPUTATION  OF  FREE  FLIGHT  AIRCRAFT  STRUCTURE  ELASTIC  DEFORMATION, 
FREE  STRUCTURE  INFLUENCE  COEFFICIENT 


The  attached  exterior  of  the  aircraft  rear  disturbance  of  balanced 

state  designed  for  free  flight  that  emerges  on  element  .1  is  f  .  ;  at 

2 .1 

this  moment  each  element  on  the  aircraft  corresponding  ineruia  force 
is  equal  to: 


*-) 


JC-  1,  2,  ••• N 


(16) 


Where  m  is  aircraft  mass;  'Iyy  is  the  pitching  rotation  inertia  of 
the  aircraft;  Xj  is  the  distance  from  the  centroid  of  element  j  to 
the  aircraft  center  mass  edge  OX  direction;  m^  is  the  mass  of  the  K 
element . 


Exterior  force  fz^  and  the  corresponding  inertia  force  system 
fT„  (K=1,2,***N)  constitute  a  balanced  force  system.  Under  the  effect 
of  these  balanced  force  systems,  displacement  of  the  i  element  of 
the  aircraft  to  the  structure  axis  system  is: 


N 

bZt“ClJz,  +  C<nflK 


(17) 


Where  c^ 
structure 


is  the  structure  influence  coefficient  of  the  corresponding 
axis  system.  With  (16)  substituting  (17)  it  makes: 


— sr 


c,„  w*— 


(18) 


to  obtain: 


(19) 
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Because  in  turbulent  motion  the  structure  axis  system  is  not  like 
the  average  body  axix  system  that  kind  of  completely  with  the  aircraft 
together  making  rigid  body  transitional  and  rotational  motion,  the 
displacement  of  each  element  in  the  aircraft,  except  for  elastic  dis¬ 
placement,  by  structure  axis  system  computations  then  cannot  avoid 
including  the  composition  of  rigid  body  displacement,  must  deduct  this 
kind  of  rigid  body  displacement,  and  then  can  obtain  the  elastic  dis¬ 
placement  of  each  element  of  the  aircraft.  Thus 

(20) 

Where  ftZ,  is  the  displacement  of  the  aircraft  related  structure  axis 
system  to  make  rigid  body  translational  motion.  Also  it  is  simply  the 
displacement  of  the  primary  point  of  the  average  body  axis  system  to¬ 
ward  the  structure  axis  system;  &£««  is  the  displacement  of  the  aircraft 
(also  is  simply  the  average  body  axis  system)  each  element  when  making 
rigid  body  rotational  motion  toward  the  structure  axis  system. 


Consider  the  condition  of  the  deformation  of  the  rear  disturbance 
edges  OX  and  0Y  axis  can  be  omitted  and  not  calculate  the  aircraft 
symmetry  of  the  left  and  right  sides.  Formulas  (14)  and  (15)  can  be 
simplified  as: 


0 


0 


Where  (20)  substitutes  (21)  and  (22)  to  obtain: 


(21) 

(22) 


With  (23)  substituting  (20)  to  obtain: 

bZt,—bZ,  m~ m*bZK — m*X**>z* 
«'  ™  1 ,  2 ,  —iV 


(23) 


(24) 
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This  is  simply  the  correct  computation  formula  for  elastic  displace¬ 
ment  of  each  element  of  the  aircraft  rear  disturbance.  The  preceding 
computation  process  clearly  shows,  by  the  displacement  of  the  related 
structure  axis  system  it  deducts  the  corresponding  rigid  body  dis¬ 
placement  &Z»,  .  Also  then  corresponding  to  the  displacement  con¬ 

version  of  the  related  structure  axis  system  is  the  displacement  of 
the  average  body  axis  system.  This  is  simply  to  say,  when  computing 
aircraft  elastic  displacement,  we  must  use  the  average  body  axis  system. 

Where  formula  (19)  substitutes  formula  (24),  it  makes: 


1  * 

C,, 

free 


Wk-ATic 


(25) 


and  we  obtain: 


bZ,, -*<(••/*> 


(26) 


Where  c.  .  „  represents  the  free  structure  influence  coefficient, 
xj  free  r 

it  expresses  the  effect  of  a  unit  of  exterior  force  that  the  element 
j  of  the  free  flight  aircraft  receives,  and  the  elastic  displacement 
in  element  i  is  produced. 


Consequently,  the  elastic  displacement  of  each  element  in  the 
aircraft  becomes: 


ibZtf)  •  Cc,;,,3  {/x,}  ( 27 ) 

V.  COMPUTATIONS  OF  FREE  FLIGHT  ELASTIC  AIRCRAFT  BALANCED  STATE  LIFT 
DISTRIBUTION  AND  REAR  DISTURBANCE  PNEUMATIC  DERIVATIVES 

Elastic  aircraft,  when  a  specified  altitude  and  specified  number 
M  make  a  level  speed  rectilinear  flight,  the  following  equation  is 
established : 


10 


(28) 


L,-i ng~  0 

A/.-  0 

{ bZi ,)  -  tat**)  {Ft,— fit  iff) 

Where  Lg  is  elastic  aircraft  lift;  Mg  is  elastic  aircraft  pitching 
motion;  F„ .  is  the  lift  of  the  effect  on  elastic  aircraft  small  ele- 
ments . 


By  equation  (28)  we  can  determine  the  angle  of  attack  <*,.»  the 
rudder  incline  angle  5  ^ ,  and  the  lift  distribution  of  the  aircraft 

this  moment  and  lift  coefficient  Cg^ g  to  deduce  resistance  coefficient 
n 

DilE  when  specified  number  M  and  specified  altitude  make  level  speed 
rectilinear  flight. 

Equation  (28)  can  be  evaluated  by  use  of  the  alternate  method. 

As  to  the  calculation  of  elastic  aircraft  rear  disturbance  pneu¬ 
matic  derivatives,  we  can  use  the  similar  method  of  computing  rigid 
aircraft  pneumatic  derivatives.  What  is  different  is  that  is  the 
computation  process  we  must  according  to  formula  (27)  seek  out  the 
attached  elastic  displacement  of  each  element  of  the  aircraft  rear 
disturbance,  and  later  again  seek  the  attached  lift  of  the  aircraft 
rear  deformation  and  pitching  motion  to  deduce  the  resistance  coeffic¬ 
ient  . 


When  calculating  the  elastic  displacement  of  each  element  of  the 
aircraft  rear  disturbance  we  must  also  use  the  alternate  method. 

In  accordance  with  the  preceding  method,  when  computing  the 
quasi-static  elastic  aircraft  pneumatic  derivatives,  because  of  ap¬ 
plied  free  structure  influence  coefficients,  we  have  already  considered 
the  use  of  exterior  force  of  disturbance  lift  and  other  corresponding 
inertia  force  systems.  Therefore  in  the  pneumatic  equation  for  air¬ 
craft  disturbance,  we  need  not  again  have  other  considerations  due 
to  the  attached  pneumatic  force  of  inertia  force  toward  the  deformation 
influence  that  is  rising.  (  Some  literature  says  that  this  section 
derivative  is  inertial  derivative.  See  the  literature  and  ^  1 0^  )  . 
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Therefore  the  corresponding  equation  for  quasi-static  elastic  aircraft 
motion  of  these  pneumatic  derivatives,  and  the  second  formula  that  is 
similar  to  the  equation  for  rigid  aircraft  pneumatic  derivative  motion 
are  computed.  We  only  must  use  the  quasi-static  elastic  aircraft 
pneumatic  derivatives  C^gFand  Cmag  with  each  rigid  body  aircraft  pne¬ 
umatic  derivatives  CT  and  C  then  can  be  replaced. 

La  ma 

VI.  ANALYSIS  OF  COMPUTED  EXAMPLES  AND  COMPUTED  RESULTS 

Using  the  preceding  methods,  we  computed  the  pneumatic  derivatives 
of  lift  distribution  of  balanced  state  and  rear  disturbance  of  cer¬ 
tain  large  scale  aircraft  on  a  TQ-6  computer  (  The  above-mentioned 
aircraft  contours  are  shown  in  the  simplified  diagram  in  Fig.  2). 


FIG.  2:  AIRPLANE  CONFIGURATION 

Fig.  3  gives  the  rigid  aircraft  and  elastic  aircraft  when  H=6000 
meters  and  M^  =0.6  under  alike  lift  coefficient  (Cg^=0.26)  when  mak¬ 
ing  level  speed  rectilinear  flight  the  wings  open  toward  the  load 
distribution  and  the  chord  toward  the  lift  distribution.  The  balance 
state  angle  of  attack  and  the  rudder  incline  angle,  elastic  aircraft 
is  °n  =3.5°,  and  81E=2.85';  rigid  aircraft  is  a1=2.76°,  and  ft,  =2.2*. 

By  this  we  can  see,  under  the  same  kind  of  conditions  of  Number  M  and 
altitude  and  flight  weight,  making  level  speed  rectilinear  flight, 
the  angle  of  attack  and  rudder  incline  angle  of  elastic  aircraft  all 
are  larger  than  those  of  rigid  aircraft. 

Fig.  4  gives  the  calculated  results  of  rigid  aircraft  and  elastic 
aircraft  when  the  preceding  altitude  and  Number  M  under  the  conditions 
of  alike  angle  of  attack  (<*=3.5* ).  We  can  see  by  the  figure,  after 
computing  the  elastic  deformation  influence,  the  aircraft  lift  coeffic- 
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ient  descends  very  much.  In  this  calculation  C^£  descends  more  than 
CL1  by  13,955  ^CL1S=0.26,  Cl1  =0.302) . 


FIG.  3:  COMPARISON  OF  LOAD  FIG.  4:  COMPARISON  OF  LOAD 
DISTRIBUTION  AND  LIFT  DIST-  DISTRIBUTION  AND  LIFT  DIST¬ 
RIBUTION  ON  RIGID  AND  ELAS-  RIBUTION  ON  RIGID  AND  ELAS¬ 
TIC  AIRCRAFT  AT  Mm  =0.6,  TIC  AIRCRAFT  AT  M  =0.6, 
Cy=0.26  ®  •  =3.5*  ® 

KEY:  (a)  Elastic  Aircraft  KEY:  (a)  Elastic  Aircraft 
(b)  Rigid  Aircraft  (b)  Rigid  Aircraft 
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Fig.  5  and  Fig.  6  separate  under  the  preceding  conditions  the 
elastic  displacement  of  the  wings  and  fuselage. 


a- 


TtSTK#) 


o-34-y 

9.2  “OT*  0.*  St*  1.0 
FIG.*’ 5:'  ELASTIC  DISPLACEMENT 


FIG.  6:  ELASTIC  DISPLACEMENT 


OF  WING.  M  =0.6,  Cy-,=0.26  OF  FUSELAGE.  M  =0.6,  Cy„=0.26 
KEY:  (a)  m&er.  S  KEY:  (a)  meter®  * 


Table  1  gives  the  computed  value  of  the  rear  disturbance  pneumatic 
derivatives  of  the  preceding  aircraft  in  the  preceding  flight  cond¬ 
itions  (H=6000meters ,  M  =o.6,  CTi=0.26).  Computed  value  (1)  is  the 
computed  result  of  rigid  aircraft.  Computed  value  (2)  is  the  comput¬ 
ed  result  of  elastic  aircraft.  By  the  table  we  can  see,  after  com¬ 
puting  the  structure  elastic  deformation  influence,  the  above-mentioned 
aircraft's  pneumatic  derivatives  have  comparatively  great  change. 


Test  Value 

! _ ' _ Computed  Value _ _ _ 

(3> 

(3)  ‘ 

Cl. 

4.05 

4.100 

4.033 

C-.  * 

-  b.osoo 

-0.500 

-0.7270 

-  0.0471 

c»„ 

0.1145 

0.117 

0.104 

.0.103 

Cl. 

0.307 

0.330  * 

0.341 

C-. 

-0.0413 

-1.0* 

-0.703 

-0.733 

CL. 

1.040 

0.7101 

0.4471 

< 

-0.310 

■  -0.110 

-7.4*1 

-7.1*0 

c». 

0.004 

o.osor 

e.osu 

CiM 

0.130 

-0.00374 

-0.0000* 

C»M 

-0.0333 

-0.01*3 

-0.0100* 

Cotm 

0.00314 

-0.00003  . 

-0.00*03 

TABLE  '1 :  M=0.6,  H=6000meters ,  CL  =  0.26,  q  =  1212  kg/m' 

U 
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For  easy  comparison.  Table  1  lists  directly  the  elastic  aircraft 
pneumatic  derivatives  that  use  the  structure  influence  coefficient 
c^j  to  calculate  deformation.  At  this  moment  the  aircraft  consolidates 
the  center  of  mass  when  there  is  no  deformation.  "Elastic"  displace¬ 
ment  is  the  structure  axis  system  of  the  aircraft  center  of  mass  that 
is  relative  to  the  primary  point  position  when  there  is  no  deformation. 
This  kind  of  method  of  computing  elastic  aircraft  pneumatic  derivatives 
in  the  past  had  been  some  means  that  were  adopted,  but  from  the  basic 
principle  we  came  to  see,  clearly  it  is  not  correct.  By  Table  1  we 
can  see  that  the  elastic  aircraft  pneumatic  derivatives  that  use  this 
method  of  computation  is  much  more  serious  than  the  derivative  that 
uses  c^  free  calculate  its  elastic  deformation.  But  this  in 
reality  is  only  a  false  appearance,  because  this  method  of  treatment 
does  not  correctly  reflect  the  realistic  condition  of  free  flight 
quasi-static  elastic  aircraft  rear  disturbance. 


VII.  SOME  PROBLEMS  OF  QUASI-STATIC  ELASTIC  AIRCRAFT  STABILITY  AND 
CONTROL 


With  the  quasi-static  elastic  aircraft  pneumatic  derivatives  that 
the  preceding  section  obtained,  .substitute  the  related  formula  of 
rigid  aircraft  stability  and  control.  We  then  obtain  the  corresponding 
formula  for  quasi-static  elastic  aircraft.  For  example,  as  to  quasi¬ 
static  elastic  aircraft  specified  speed  focus,  we  have: 


-  - 

igf- 


As  to  specified  load  focus,  we  have: 


(29) 


(30) 


As  to  specified  normal  level  rectilinear  flight,  the  rudder  incline 
angle  with  the  variable  frequency  of  the  the  velocity,  we  have: 


(31) 
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Where  I—,  and  ZwiUv  the  X  coordinates  of  the  aircraft  centroid, 
the  specified  speed  focus,  and  the  specified  load  focus,  and  their 
ratio  to  the  average  pneumatic  chord.  (A&/A9)a  is  the  ratio  of  the  in¬ 
crements  of  the  rudder  incline  angle  in  level  speed  rectilinear  flight 
and  the  increments  of  flight  velocity,  in  which  A*"  »/»«•.. 

In  Table  1,  the  value  of  the  quasi-static  elastic  aircraft 
corresponding  pneumatic  derivatives  substitute  formulas  (29),  (30), 
and  ( 31 )  to  obtain: 

-0.177,  0.1416,  (A6/AS)a-0.1 

As  to  rigid  aircraft,  the  preceding  coefficient  values  are: 

-0.1892,  X.f-xJ-,- -0.18,  A&/A5-0.095. 

Thus  we  can  see,  after  computing  the  elastic,  deformation  influence 
the  specified  speed  focus  '  and  the  specified  load  focus  both  move 
forward,  and  the  value  of  A&/A9  yet  has  increases. 

In  order  to  understand  the  pneumatic  characteristics  of  quasi¬ 
static  aircraft  high  subsonic  speed  pressure  time,  we  still  computed 
the  focal  position  of  elastic  aircraft  and  rigid  aircraft  when  M  =0.8; 
its  value  is  as  Table  2  expresses: 


, 

Rigid  Aircrai 

r 

Elastic 

Aircraft 

24(8 

21(5 

24(1 

X.,  -  X.. 

-0.233 

-0.202 

-  0.199 

-0.22 

-  0.1917 

I 

-  0. 188 

TABLE  2 


The  preceding  computations  explain  that  when  M  increases  from 
0.6  to  0.8,  and  speed  pressure  q  increases  from  121.2  kg/m  to  2155  kg/m 
then  elastic  aircraft  specified  speed  focus  moves  back  in  volume  less 
than  rigid  aircraft  does.  This  is  because  the  influence  of  the  speed 
force  rising  elastic  deformation  lowers  the  influence  of  a  section  of 


m 


Number  M.  By  Table  2  we  can  yet  see  that  when  Number  M  is  fixed, 
with  the  speed  pressure  increase,  the  specified  speed  focus  and  the 
specified  load  focus  both  move  forward,  and  specified  speed  and  spec¬ 
ified  load  degree  of  stability  both  descend. 

Because  the  factors  that  effect  elastic  aircraft  pneumatic  char¬ 
acteristics  are  comparatively  complex,  by  certain  opinions  or  con¬ 
clusions  of  aircraft  computed  results  that  have  been  gained,  no  cer¬ 
tainties  can  be  extended  to  other  aircraft.  For  example,  in- literature 
(6)  the  value  of  the  specified  load  focus  of  a  certain  aircraft  when 
M_=0.8  is  claculated.  The  above  text  considers  that  when  M  =0.8, 

CD  CD 

with  speed  pressure  increases,  specified  load  focus  position  contin¬ 
ually  moves  back,  when  speed  pressure  is  very  large,  the  specified 
load  focus  is  far  behind  the  aircraft  centroid.  Moreover,  the  var¬ 
iations  are  not  sensitive  to  the  centroid.  But  in  this  text,  the 
conditions  of  the  preceding  computations  are  not  like  this.  When 
M^O.8,  with  speed  pressure  increase,  specified  load  focus  on  the 
contrary  moves  forward.  Therefore  we  consider  that  if  the  literature 
f6)#  computed  results  are  believable,  they  only  represent  particular 
conditions  of  certain  types  of  aircraft.  Moreover,  we  cannot  obtain 
fromthis  a  general  conclusion. 

Lastly  we  want  to  briefly  discuss  some  problems  that  merit  at¬ 
tention  when  using  some  other  methods  to  treat  quasi-static  elastic 
aircraft  dynamic  state  characteristics.  In  recent  years,  foreign 
documents  have  brought  out  some  methods  to  use  the  so-called  "zero 
quality  pneumatic  derivatives"  and  "inertia  derivatives"  to  compute 
elastic  aircraft  pneumatic  characteristics.  The  so-called  "zero 
quality  pneumatic  derivative"  is  simply  on  the  aircraft  rear  distur¬ 
bance,  and  only  computes  the  elastic  aircraft  pneumatic  derivatives 
of  the  exterior  force  rising  "structure  deformation"  gains.  The  so- 
called  "inertia  derivative"  is  simply  on  the  aircraft  rear  disturbance, 
and  only  considers  the  effect  of  the  derivative  of  the  aircraft  inertia 
force  rising  "structure  deformation"  gains.  When  taking  the  exterior 
force  system  and  its  corresponding  inertial  force  system  to  separately 

•  Computed  formulas  in  literature  6  of  specified  load  focus  are  not 
correct.  See  the  following  discussion  in  this  text  or  literature  4  . 
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compute  the  structure  deformation  that  has  arised,  because  they  both 
are  not  balanced  force  systems,  when  computing  deformation,  the  sus¬ 
taining  reaction  force  that  is  used  is  not  equal  to  zero.  Therefore, 
in  "zero  quality  pneumatic  derivatives"  and  "inertia  derivatives" 
computations,  the  sustaining  effect  exists.  When  the  sustaining  con¬ 
ditions  change,  the  numeric  value  of  the  derivative  also  changes. 
Therefore,  by  themselves  "zero  quality  pneumatic  derivatives"  and 
"inertia  derivatives"  are  not  of  realistic  significance.  But  after 
substituting  elastic  aircraft  pneumatic  equations  with  these . deriva¬ 
tives  (This  kind  of  quasi-static  elastic  aircraft  pneumatic  equations 
and  rigid  aircraft  pneumatic  equations  are  not  the  same),  because 
all  the  exterior  force  systems  of  disturbance  rising  and  their  cor¬ 
responding  inertia  force  still  constitute  balanced  force  systems,  the 
sustaining  effect  of  these  derivatives  mutually  reduce.  Using  this 
kind  of  equation  and  derivative  to  compute  the  dynamic  state  charac¬ 
teristics  of  quasi-static  .elastic  aircraft,  its  conclusion  is  still 
correct.  But  the  "zero  quality  pneumatic  derivative"  alone,  due  to 
the  sustaining  influence  included  among  it,  cannot  be  used  to  compute 
and  analyze  the  control  and  stability  characteristics  of  the  air¬ 
craft.  Therefore,  any  other  derivative  formula  using  "zero  quality 
pneumatic  derivatives"  or  uses  when  computing  aircraft  elastic  de¬ 
formation  not  balanced  force  system  gains  to  research  the  fixed  speed 
focus  of  elastic  aircraft  and  the  fixed  load  focus  and  its  control  and 
stability  performance  generally  speaking  is  not  correct.  Literature 
<^6)when  computing  the  fixed  load  point  of  elastic  aircraft,  the  pneu¬ 
matic  derivative  that  is  used  clearly  brings  on  the  sustaining  influence, 
because  what  it  uses  when  computing  aircraft  structure  deformation  is 
not  balanced  force  system.  Therefore  the  sustained  reaction  force  is 
not  equal  to  zero.  At  the  same  time,  the  coordinate  system  that  the 
above  text  uses  when  computing  the  aircraft  structure  deformation  is 
not  a  balanced  axis  system.  This  kind  of  computed  results  are  not 
accurate . 

In  some  other  literature  (for  example  7  and  10)  some  similar  prob¬ 
lems  also  exist.  Because  literature  4  has  already  made  detailed 
discussion  of  these  problems,  we  won't  repeat  it  here. 
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This  text  has  already  been  checked  and  approved  by  Prof.  Zhang 
Guilian,  who  also  pointed  out  valuable  opinions.  We  hereby  express 
our  thanks. 
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PRELIMINARY  STUDY  ON  VARIABLE  POROSITY  WALLS  FOR  A  TRANSONIC 
WIND  TUNNEL 


Variable  Porosity  Wall  Research  Group,  Nanjing  Aeronautical  Institute 
Penned  by  Zhang  Qiwei 

Abstract : 

In  order  to  reduce  the  wall  interference  and  improve  quality  of  the  flow 
field  in  a  transonic  wind  tunnel,  a  set  of  variable  porosity  walls  with  60‘  incli¬ 
ned  holes  have  been  designed  and  manufactured.  The  open-area  ratio  of  the 
walls  can  vary  continually  from  zero  to  9.2  percent.  The  walls  have  been  used 
in  a  600mm X  600mm  tran-  and  supersonic  wind  tunnel  with  solid  side  walls. This 
paper  describes  the  general  characteristics  of  the  variable  porosity  walls  and  the 
preliminary  results  of  calibration  at  Mach  numbers  ranging  from  0.6  to  1.2. 


SYMBOLS: 

M  Airflow  Mach  number 

Mt  Test  section  nucleus  uniform  internal  area  average  Mach  number 

Mo  Computed  airflow  Mach  number  of  static  pressure  from  station 
reference  point 

Mw  Computed  airflow  Mach  number  of  static  pressure  from  test 
section  side  wall 

AM.  AM.-M.-M. 

AM.  AM.-(M--M.)with  model  — (M»— M.)  without  model 
Re  Unit  length  Reynold's  number  (1/m) 

P  Static  pressure 

Po  Front  chamber  total  pressure 

®  Wall  plate  expansive  angle  (expanse  is  straight) 

K  Open  area  ratio  of  porosity  wall 

d  Diameter  of  porosity  wall 

A  Dynamic  plate  shift  distance  of  variable  porosity  wall  plate 
D  Diameter  of  20*  cone-cylinder  model  cylinder  section 

H  Wind  tunnel  test  section  altitude 

La  Airflow  acceleration  section  length 

x  Distance  from  test  section  entrance  cross  section 

x.j  Distance  from  top  point  of  20°  cone-cylinder  model 
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ffw  Mean  root  error  of  the  test  section  nucleus  flow  uniform 
internal  area  Mach  number 

Cy  Derivative  (l/o  )  of  the  zero  facing  angle  time  lift  cooefficient 
to  the  facing  angle 

xp  Ratio  of  the  the  distance  of  the  zero  facing  angle  time  focus 
from  the  top  point  of  the  model  and  the  largest  diameter  of 
the  model  fuselage 

CXq  Resistive  coefficient  of  the  zero  facing  angle  time 
£Cxq  Resistive  coefficient  increments,  £CXq=CXq-Cxq (M_q  g) 

I .  PREFACE 

In  order  to  fully  weaken  the  interference  of  the  tunnel  walls 
of  the  transonic  wind  tunnel,  outside  the  country  many  wind  tunnels 
use  variable  porosity  walls ^  .  in  order  to  improve  the  flow  field 

characteristics  of  the  transonic  wind  tunnel  and  understand  some 
particular  problems  of  using  variable  porosity  walls,  we  have  designed 
and  manufactured  on  our  own  a  set  of  continuable  variable  porosity 
walls.  In  1980  this  set  of  porosity  walls  was  used  in  a  tunnel  with  a 
high  speed  of  600mmx600mm  to  conduct  preliminary  flow  field  calibra¬ 
tion:  researched  in  the  scope  of  M=0. 6~1 . 2  porosity  and  the  influence 
of  wall  plate  expanse  angle  to  the  flow  field,  and  calibrated  the 
interrelation  between  the  station  reference  point  Mach  number  and  the 
nucleus  flow  uniform  internal  area  average  Mach  number,  using  the  20* 
porosity  cone-cylinder  model  to  conduct  calibration  pressure  tests. 

It  was  preliminarily  determined  to  use  from  new  on  the  scope  of  &  and 
K,  and  also  conducted  AGARD-B  standard  model  calibration  pressure  tests. 
This  text  explains  the  general  situation  of  this  kind  of  variable  por¬ 
osity  wall  and  preliminary  calibration  results. 

II.  SIMPLIFIED  CONDITIONS  OF  A  WIND  TUNNEL 

This  wind  tunnel  is  a  direct  flow  down  temporary  formula  wind 
tunnel.  Its  minute  details  are  seen  in  material  {^3}.  The  test  section 
model  cross  section  is  a  square  of  600mmx600mm.  The  test  section  length 
is  1  580mm.  The  test  section  Mach  number  range  is  0.5vO.5.  When  M  =  0.5"^ 
-'1.2  a  solid  speed  of  sound  nozzle  is  used,  depending  on  changing  the 
front  chamber  total  pressure  to  change  the  Mach  number.  The  left  and 


21 


T 


*7* 


right  side  walls  of  the  test  section  are  fixed  parallel  solid  walls, 
which  are  porous  from  top  to  bottom,  and  the  expanse  angles  can  be 
adjusted.  The  station  depth  is  comparatively  small.  The  most  shallow 
place  is  only  170mm.  Moreover,  there  is  no  auxiliary  air  exhaust 
system.  As  our  first  step  we  only  take  the  straight  porosity  walls 
that  were  formerly  used  and  switch  over  to  variable  porosity  walls 
to  conduct  tests.  When  M=1. 5^3.5  this  wind  tunnel  uses  replaceable 
two  dimensional  nozzles  to  change  the  Mach  number,  and  the  test  sec¬ 
tion  walls  are  all  solid  walls. 

III.  GENERAL  CONDITIONS  OF  VARIABLE  POROSITY  WALLS 

Variable  porosity  walls  are  composed  of  steel  plates  that  have 
two  layers  open  with  holes.  Detailed  structure  can  be  seen  in  mater¬ 
ial  £ 4),  and  a  simplified  diagram  is  shown  in  Fig.  1.  The  layer  of 
porosity  plate  near  the  airflow  is  fixed  and  doesn't  move.  The  other 
layer  of  porosity  plate  can  shift  back  and  forth.  When  the  holes  of 
the  two  plates  are  completely  aligned,  the  wall  plate  porosity  gets 

the  largest.  According  to  the  computations  in  Fig.  1  ,  K.,  .  = 

2  lar  g  ©  s  x> 

=  (jtx4.5  )/ ( 65x1  0 )  =  9 . 8% ,  but  because  the  fixed  position  screw  occupies 
a  separate  hole  position,  the  realistic  largest  porosity  is  9.2%. 


FIG.  1  DETAILS  OF  VARIABLE  POROSITY  WALLS 

(a)  Simplified  diagram  of  variable  porosity  wall 

(b)  Variable  porosity  test  section 

KEi:  (1)  Airflow;  (2)  Fixed  plate;  (3)  mobile  plate 
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When  the  holes  of  the  two  plates  are  completely  spread  apart,  the 
porosity  tends  toward  zero.  The  holes  are  evenly  distributed,  the 
back  of  the  moving  plate  of  the  front  wall  can  plug  into  the  block¬ 
ing  plate  of  various  models,  to  change  the  airflow  acceleration  char¬ 
acteristics.  Experiments  demonstrate  that  when  not  adding  any  block 
plate  the  airflow  acceleration  section  is  shortest  (when  K=6?,  La/H 
can  shorten  to  0.4),  and  at  this  moment  in  the  uniform  area  the  air¬ 
flow  fluctuation  also  is  not  great.  Therefore,  now  we  do  not  add  a 
blocking  plate. 


The  model  support  of  this  wind  tunnel  is  located  in  the  inside 
of  the  test  section.  When  wall  plate  porosity  varies  very  little, 
the  transonic  support  plug  up  effect  will  be  extremely  severe.  There¬ 
fore  on  the  support  area  wall  plate  is  added  ten  cracks  that  are  5mm 
wide  and  200mm  long  (see  the  photograph  in  Fig.  1). 


The  moving  plate  of  this  wall  plate  is  run  by  a  model  DG-6  elec¬ 
tronic  motor.  When  reducing  porosity,  the  moving  plate  contrary  air 
flow  direction  shifts.  Its  shift  distance  is  measured  by  model  LVDT 
displacement  sensor.  Measurement  precision  is  +0.1mm.  The  relation¬ 
ship  between  the  moving  plate  shift  distance  A  and  porosity  K  is  ac¬ 
cording  to  the  following  computations  (when  we  make  K=K^a  A=0 ) : 


K~K**X  Hrf"«  "kirH-a-y 1 

*  « 


IV.  CONTENTS  OF  EXPERIMENTS 

n 

Emperiment  Mach  number  scope  is  0.6--1.2,  Re=(1.3  2.3)x10  /m; 
Porosity  variation  scope  is  0— *7^ ;  Wall  plate  expanse  angle  variation 
scope  is  0-~0.6  .  Front  chamber  total  pressure,  station  pressure,  and 
test  section  sideplate  static  pressure  is  calibrated.  We  use  fixed 
in  the  test  section  on  the  center  line  axial  static  pressure  sounding 
tube  to  calibrate  test  section  nucleus  flow  Mach  number  distribution. 
The  axial  static  pressure  sounding  tube  external  diameter  is  28mm; 
on  the  tube  each  interval  of  40mm  has  a  calibration  pressure  hole, 
and  top  to  bottom  staggered  arrangement.  We  also  calibrated  the  sur- 
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face  distributed  pressure  of  a  20°  cone-cylinder  model.  The  length 
of  the  cone-cylinder  model  is  720mm;  the  diameter  of  the  cylinder 
section  is  68mm;  the  plug  proportion  weight  is  1/5.  The  position  of 
the  axial  sounding  tube  and  the  20°  cone-cylinder  model  in  the  test 
section  is  shown  in  Fig.  2.  Aside  from  this,  we  still  conducted 
AGARD-B  standard  tests  to  measure  pressure  under  different  porosity 
variation  regulations. 


FIG.  2  LOCATION  OF  CENTERLINE 
STATIC  PRESSURE  AXIAL  SOUNDING 
TUBE  AND  20°  CONE-CYLINDER  MOD¬ 
EL  IN  THE  TEST  SECTION 
KEY:  (a)  Calibration  hole  area 

V.  PRINCIPAL  TEST  RESULTS 

1 .  The  Flow  Field  of  Nucleus  Flow 

Figures  3  and  4  separately  give  the  Mach  number  distribution 
of  the  nucleus  flow  under  typical  conditions  and  the  (7^  value  of  this 
wall  plate  under  general  conditions. 

(1)  The  influence  of  &  to  the  nucleus  flow  field:  when  &  value 
increases,  subsonic  speed  airflow  uniform  degree  gets  better,  super¬ 
sonic  speed  airflow  acceleration  gets  faster.  If  b  enlarges,  then 
excessive  acceleration  emerges  on  the  supersonic  speed  area,  and  causes 
the  flow  field  uniform  degree  to  get  worse,  b  decreases  then  it  causes 
the  Mach  number  to  be  unable  to  reach  1.2.  6  can  be  chosen  and  used 
in  the  range  of  0.2*^  0.5*. 
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(a)  Supporting  area  without  (b)  Supporting  area  with  crack 
crack*  =0.2*  *=0.5* 


FIG.  3  TYPICAL  CENTERLINE  MACH  NUMBER  DISTRIBUTION 


FIG.  4  CENTERLINE  LOCAL  MACH  NUMBER 
DEVIATION  (RMS) 

KEY:  (a)  Symbol;  (b)  Straight  hole  wall 


(2)  The  influence  of  K  to  the  nucleus  flow  field:  When  the  value 
of  K  increases,  subsonic  speed  airflow  uniform  degree  gets  better,  and 
supersonic  speed  airflow  acceleration  gets  faster.  If  K  enlarges,  then 
excessive  acceleration  emerges  in  the  supersonic  speed  area,  and  causes 
airflow  uniform  degree  to  get  worse.  From  the  flow  field  uniform 
degree  we  come  to  see  that  when  subsonic  speed  takes  K=6^7%  it  is 
better,  at  supersonic  speed  the  K  value  is  best  to  change  with  the 
Mach  number,  its  numeric  value  is  related  to  &  ,  when  b  is  larger 


25 


'  *.'**«-  ■  -  ;  - 


we  should  apply  a  smaller  K  value.  For  example,  when  Mt=1.20, 
we  takefc=0.5a,  and  K=6$  is  better;  when  Mt^1.15,  if  5=0. 5°  »  then 
to  take  5$  is  better;  if  b-=0.2  ,  then  to  take  K$5'*6%  is  better. 

(3)  Influence  of  supporting  region  wall  plate  slot:  The  rigid 
support  of  this  wind  tunnel  is  located  in  the  test  section.  Its  plug 
ratio  weight  is  4.26$.  By  Fig.  3  we  can  see,  if  the  supporting  region 
wall  plate  has  no  slot,  then  when  porosity  is  small  the  influence  of 
the  supporting  plug  effect  to  the  upper  area  is  very  serious,  causing 
the  airflow  uniform  region  to  greatly  shorten,  and  the  above-mentioned 
phenomenon  behind  the  supporting  region  slot  greatly  lightenes.  We 
can  see  that  the  supporting  region  wall  plate  slot  is  neccessary. 

To  sum  up,  the  flow  field  uniform  degree  of  this  wall  plate  is  better 
than  the  fixed  straight  hole  wall  that  this  wind  tunnel  formerly  used 
(Fig.  4). 


2.  The  Station  Reference  Point  Mach  Number  Me 

Through  calibration  we  discovered  that  the  pressure  distribution 
of  the  lower  part  of  the  station  is  relaticely  uniform,  therefore  we 
each  took  a  calibration  pressure  point  of  comparatively  advanced  pos¬ 
ition  in  the  lower  section  up  and  down  the  station  and  according  to 
the  total  made  a  reference  point.  The  same  time  as  calibrating  the 
nucleus  flow  Mach  number,  we  calibrated  the  reference  point  Mach  num¬ 
ber  Me.  Figure  5  gives  the  deviation  value  of  Me  and  Mt  under  typ¬ 
ical  conditions. 

When  the  supporting  region  wall  plate  has  a  slot,  the  effect  of 
the  variation  toward^Mc  is  relatively  small.  In  Fig.  5  (a)  the  scope 
of  variation  of  Me  value  approximately  corresponds  to  2°"^,  this  and 
the  ADP  wall  plate  test  results  in  literature  (fn  are  similar.  When 
the  wall  plate  has  no  open  slot,  with  the  reductio  of  the  K  value, 

AMc  markedly  reduces,  this  and  the  LAR  wall  plate  test  results  of 
literature  ( 5}  are  similar.  We  can  see  that  the  wall  plate  slot 
(slot  ventilation  area  is  not  following  K  variation)  can  reduce  the 
effect  of  K  on  AMc. 

In  order  to  research  the  effect  of  the  model  existence  on  M,  when 
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FIG.  5:  TYPICAL  /VMc  VALUES 


FIG.  6:  MODEL  EFFECT  ON  THE  CALIBRATING 
MACH  NUMBERS  (8  =  0.5*,  supporting  region 
with  slot) 

conducting  20  cone-cylinder  model  calibration  tests,  we  also  cal¬ 
ibrated  the  model  front  of  Mw  (when  M^l.05,  the  calibration  point 
is  located  in  the  center  of  the  side  wall  of  155mm  of  the  model  top 
point  front,  when  M,>1.05  Mw  takes  the  side  wall  calibrated  Mach  num- 
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ber  of  the  front  of  the  model  top  shock  wave).  1'hrough  computing 
Md=(Mw-Mc)w.th  fflodel-(Mv-Mc)w.thout  model  of  the  corresponding 
calibration  point  that  is  obtained,  such  as  Fig.  6  expresses. 

As  to  the  20*  cone-cylinder  model  where  the  plug  ratio  is  1%, 
notice  that  A^d  is  not  very  large. 

3.  20 Cone-Cylinder  Model  Calibration  Test  Results 

General  test  results  are  shown  in  Fig.  7.  In  the  figure,  the 
solid  line  is  the  non-interence  curve  (refer  to  the  photograph  in 
material  (6));  The  test  point  corresponding  M  is  the  Mt  that  is  arriv 
ed  at  by  Me  seen  in  Fig.  5;  not  carried  out  AMd  correctly.  Because 
the  tail  supporting  shaft. of  this  model  is  too  short,  in  the  figure 
the  position  of  x.j/D=7  is  equal  to  x=900mm,  making  the  model  back 
section  of  x^/E >7  receiving  at  different  degrees  the  supporting  ef¬ 
fect  (Refer  to  Fig.  3).  Hereafter  we  will  draft  a  longer  tail 
supporting  shaft  to  conduct  tests. 


FIG.  7:  PRESSURE  DISTRIBUTION  ON  THE  20®  CONE-CYLINDER  MODEL 
KEY:  (a)  symbol 
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When  M<1 ,  K  =  5  7%  both  can  obtain  relatively  good  pressure  distribu¬ 
tion  that  is  with  a  non-interference  curve  (in  agreement  with  lit¬ 
erature  (l^  and  £2^). 

When  Mkl ,  the  expansion  wave  tunnel  wall  reflection  back  that 
the  model  shoulder  section  sends  out  becomes  a  compressed  wave  and 
reaches  on  the  model  (When  M=1  it  reaches  x^/D=5.5  approximately, 
when  M=1.15  it  reaches  x^/D=9.5  approximately).  This  shows  that  the 
porosity  has  enlarged.  With  the  reduction  of  K  value,  this  reflection 
wave  weakens.  The  results  when  M=1  and  of  Fig.  9.21  of  reference 
material  (6^1  are  close.  But  when  number  M  is  comparatively  large, 
a  four  wall  open  hole  wind  tunnel  emerges  'in  which  there  are  no 
contradictory  phenomena.  For  example  when  M=1.1,  with  the  reduction 
of  K  value , approximately  x^/D=4  emerges  obvious  compression.  This 
is  the  result  of  the  front  section  shock  wave  longitudinal  tunnel  wall 
reflexionbecoming  compression  wave  arriving  on  the  model.  This  makes 
it  clear  to  the  front  section  shock  wave,  the  K  value  is  also  too  small. 
In  addition,  when  M=1.2  and  when  K  increases  to  7%  it  still  assumes 
a  solid  wall  effect.  We  consider  the  primary  reasons  that  these  phe- 
nomina  emerge  are  probably  that  this  wind  tunnel  only  has  two  wall 
holes.  When  M=1^1.05,  the  K  requirements  to  eliminate  the  waves  are 
very  small,  the  wave  reflection  of  side  wall  production  is  not  strong. 
With  the  increase  of  the  Mach  number,  the  K  value  that  is  required  to 
eliminate  waves  is  also  very  large,  but  the  solid  wall  effect  of  the 
two  sides  then  becomes  more  and  more  striking.  In  the  tests  we  also 
calibrated  pressure  distribution  of  certain  points  on  each  line  of 
iheleft  and  right  sides  of  the  20*  cone-cylinder  model  and  discovered 
that  the  pressure  distribution  of  each  line  form  top  to  bottom  is 
basically  alike.  When  the  K  value  of  the  whole  wall  changes,  the 
pressure  distribution  of  each  line  simultaneously  changes.  From  the 
reflection  wave  system  of  the  four  walls,  before  arriving  on  the  mod¬ 
el,  we  see  that  they  already  passed  through  a  series  of  complex  mixed 
processes  (but  this  is  not  a  simple  counteraction  relationship). 

4.  AGARD-B  Standard  Calib ’ation  Test  Results 

In  order  to  make  clear  the  effect  of  variable  porosity  on 
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conventional  calibration  tests,  we  yet  carried  out  three  kinds  of 
AGARD-B  standard  calibration  tests The  plug  ratio  of  the  model 
in  the  test  section  is  0.71*,  The  wall  plate  expanse  angle  is  0,2* , 
porosity  variation  regulations  are  shown  in  Fig.  S.  Regulation  A  is 
with  porosity  fixed  at  6%,  regulation  C  is  the  proposed  applied 
porosity  value  based  on  previous  flow  field  calibration  and  20*  cone- 
cylinder  model  test  results.  Regulation  B  is  situated  between  A  and 
C. 

Typical  test  results  are  shown  in  Fig.  9.  In  order  to  make  com¬ 
parison,  Fig.  9  separately  draws  the  test  results  of  this  tunnel  in 
1974  using  K=22.4£  straight  hole  wall  and  the  combined  calibration 
results  of  each  wind  tunnel  inside  the  country  and  foreign  standard 
model  calibration  curves.  By  the  figure  we  can  see: 

(1 )  The  porosity  wall  test  results  of  K= 6%  and  the  straight  porosity 
wall  test  results  of  K=22..4?  are  balanced. 


BER  KEY:  (a)  Foreign  calibration  curve; 

(b)  China  calibtation  curve;  (c)  Vari 
able  porosity  regulation;  (d)straight 
hole  wall;  (e,f,g)  hole  wall. 


30 


(2)  When  porosity  is  close  to  regulation  C,  ^xq  is  clearly  close 
to  no  plug  interference  value.  The  reason  can  possibly  be:  when 
K  value  is  comparatively  large,  the  expanded  wave  longitudinal 
tunnel  wall  reflection  that  the  model  front  section  sends  out  is 
compressed  and  arrives  on  the  model  back  section  (See  Fig.  7),  using 
the  pressure  lift  of  the  wing  back  section,  and  resistance  reduces. 
With  the  reduction  of  K  value,  the  above-mentioned  reflective  wave 
weakens.  When  we  arrive  at  regulation  C,  the  model  back  section,  al¬ 
though  not  completely  at  a  non-interference  state,  the  model- front 
section  pressure  has  raised  (See  Fig.  7).  Both  have  a  comprehensive 
effect,  causing  resistance  to  increase  to  non-interference  value. 

In  the  calibration  test  suitable  to  variable  porosity  we  see  that  we 
can  reduce  the  tunnel  plug  interference  when  M^I . 

(3)  The  influence  of  porosity  variations  to  Cy  and  Xp  all  are  not 
noted.  The  reason  could  possibly  be  that  AGARD-B  standard  test  is 
not  sensitive  to  tunnel  wall  lift  disturbance,  and  there  is  no  tail 
wing. 

VI.  CONCLUDING  REMARKS 

We  in  the  up  and  down,  two  walls,  with  holes  transonic  wind 
tunnel  test  section  carried  out  prliminary  experimental  research  on 
variable  porosity  walls.  Through  experiments  we  preliminarily  grasped 
some  distinguishing  features  of  the  flow  field  when  using  variable 
porosity  walls.  The  test  results  express  variable  porosity  walls  are 
beneficial  to  reduce  tunnel  wall  plug  interference  when  M-^1 .  But  be¬ 
cause  there  are  only  two  walls  with  holes,  when  M^»1  it  cannot  elimin¬ 
ate  the  reflection  wave  system  of  the  right  and  left  side  walls,  and 
in  the  20“  cone-cylinder  model  calibration  pressure  test  eliminating 
the  wave  effect  is  not  ideal.  We  intend  from  now  on  to  further  dev¬ 
elop  the  research  on  wing  model  calibration  pressure  tests. 
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A  NEW  SIMULTANEOUS  ITERATION  ALGORITHM  FOR  EVALUATION  OF 
EIGENPROBLEMS  IN  LARGE  STRUCTURES 


Liu  Guoguang  and  Li  Junjie,  Aircraft  Structural  Mechanics  Research 
Institute 

Abstract : 

A  new  improved  algorithm  of  simultaneous  iteration  is  presented  for  eva¬ 
luation  of  generalized  eigenproblems  in  dynamic  analysis  of  large  structures. The. 
convergence  of  this  algorithm  is  proved  by  the  concepts  of  Ea  subspace  and 
eigendirection  and  some  details  of  how  to  perform  this  algorithm  in  computer 
are  discussed. 

Based  on  statistics  of  computational  quantity  required  for  each  step  of  itera¬ 
tion  and  a  great  number  of  practical  computations,  it  is  shown  that  in  conver¬ 
gence  rate  the  present  method  is  comparable  to  some  simultaneous  or  subspace 
iteration  algorithms  available  up  to  now,  but  it  can  cut  down  the  computer  run 
time  for  each  step  of  iteration  and  raise  the  computation  efficiency. 

I.  FORWARD 

The  eigenproblem  in  structure  dynamic  analysis,  after  passing 
through  some  strange  treatment  (for  example  primary  point  displace¬ 
ment)  its  model  is: 


Kx-XAfx 


(D 


In  which  K  and  M  are  n  exponent  solid  symmetry  positive  fixed  and 
scattered  distance.  The  eigenvalue  of  the  present  design  system  (1) 
is  A„  .  Moreover,  there  is: 

0<A,<A,< . <JU 


The  corresponding  eigen  vectors  are  x^ ,  x2» . xno* 

As  to  large  structures,  the  exponent  number  of  system  (1)  can 
be  raised  to  several  thousand.  Under  this  kind  of  condition,  using 
the  scattered  characteristics  of  K  and  M  will  be  extremely  necessary. 


T 


Therefore,  simultaneous  iteration  algorithm  already  becomes  one  of 
the  strong  powerful  methods  of  evaluating  the  section  lower  level 
eigen  face  (eigen  value  and  corresponding  eigen  vector)  of  system  (1). 
Previously,  the  structure  algorithm  method  put  forward  by  Rutishauer, 
McCormick,  Bathe,  and  Nicolai  was  prevelent.  Among  them  the  previous 
two  kinds  of  algorithm  were  conducted  b.y  Cholesky  on  resolution  to 
M,  later  again  turned  system  (1)  into  a  standard  eigenproblem  Ay=Vy 
solution.  If  we  don’t  consider  the  production  error  influence  by 
resolution  of  M,  using  this  kind  of  algorithm  solution  system  (1) 
section  eigenproblem  also  is  not  beneficial  (See  Table  4).  This  is 
not  only  because  the  resolution  of  M  in  itself  will  waste  a  certain 
amount  of  machine  time,  but  also  because  resolution  "pack"  will  be 
able  to  produce  new  non-zero  elements,  and  therefore  after  it  is  used 
each  step  of  iteration  increases  a  certain  amount  of  computer  run  time, 
thereby  effecting  computer  efficiency. 

Additionally,  the  algorithm  of  Reinsch  and  Nicolai  used  Gram- 
Schmidt  perpendicular  process,  and  therefore  we  can  stop  the  already 
weakened  test  vector  and  not  again  add  later  iteration.  This  treat¬ 
ment  technique1'  ^is  generally  called  the  "reduced  dimension"  treatment. 
When  the  number  of  treated  eigen  faces  is  comparatively  large,  using 
the  reduced  dimension  treatment  can  greatly  raise  computer  efficiency 
(See  Table  6). 

This  text  gives  a  new,  improved  algorithm  and  is  comparative  to 
the  preceding  algorithm,  aside  from  the  unresolved  M  and  the  useable 
reduced  dimension  treatment,  it  yet  has  the  primary  characteristic  of 
further  reducing  each  step  of  iteration  that  wastes  computer  run  time. 

In  the  following  section  we  will  simply  explain  the  Ek  subspace 
concept  of  system  (1)  to  prove  the  basic  theory  of  convergence  pro¬ 
vided.  The  third  section  explains  some  certain  particular  problems 
of  algorithm  and  recounts  some  aspects  realized  in  the  computer. 

Method  comparisons  and  numeric  value  test  results  will  be  given  sep¬ 
arately  in  the  two  sections  4  and  5, 


34 


We  will  use  capital  letters  to  express  matrix.  Especially,  we 

will  use  R  to  express  the  upper  triangle  matrix,  and  use  r^.  to 

express  its  (i,  j)  elements.  We  will  use  D  to  express  the  matrix  to 

the  angle,  and  d . .  to  express  its  j  diagonal  elements.  We  will  use 

X,  Y,  Z,  and  S  to  express  nxp  vector  matrix,  and  separately  use  x., 

y.,  z.,  and  s.  to  express  their  j  series  of  vectors.  We  will  use  the 
#  J  J 

capital  letter  0  to  express  the  same  exponent  small  mass. 


II.  THE  SUBSPACE  OF  Ek  IN  SYSTEM  (1) 


Assuming  x^°^,  x^^ 


as  in  the  n  dimension  subspace 


holding  a  linear,  non  related  vector  form.  If  EQ  is  the  subspace 
that  they  opened,  then  the  Ek  subspace  of  system  (1)  can  be  defined  as: 


(3) 

By  the  positive  fixed  symmetry  of  K  and  M  we  can  know  that  the 
dimension  of  the  Ek  subspace  to  the  arbitrary  k  all  is  p.  Moreover 
lira  2k  exists,  and  is  an  unvariable  subspace  of  system  (1).  If  it  is 

by  the  front  of  system  (1)  p  lower  level  eigen  vector  x^ ,  . ,  x 

opened,  then  it  calls  Ek  to  be  a  stable  convergent. 


In  this  theory,  under  stable  convergent  conditions,  on  Ek  all 

we  can  find  to  call  it  is  the  best  direction  of  p  vector  u^ . . 

Up  and  we  have : 


IK-Xfll-CKg'?) 


(4) 


Among  which  (5) 

Because  there  is  no  method  to  compute  the  best  direction,  every  al¬ 
gorithm  of  the  simultaneous  iteration  is  p  vector  on  Ek  that  goes 

through  a  form  already  known  to  compute  a  new  vector  xjk^, .  x^, 

to  use  when  K  is  amply  large,  they  gradually  tend  toward  the  favorable 
direction.  Therefore: 

||*J« ( 

Here  q^.  is  the  same  as  in  formula  (5). 
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The  preceding  simultaneous  iteration  various  algorithms,  when 


computed  on  Ek  have  the  preceding  quality  x) 
different . 


(k) 


Each  method  used  is 

The  algorithm  this  text  proposes  will  adopt  of  Ek  the 
computation  system: 


(7) 


the  eigen  vector  Xk=(x 


(k) 

1 


.(k) 


),  it  satisfies  that: 


XjAfX»=*J, 


Xf(KAf'!OX*«Di 


(8) 

(9) 


In  which  Dk  is  p  exponent  matrix  to  the  angle,  and  itself  from  start 

( k) 

to  finish  assumes  its  diagonal  element  d.  .  to  satisfy  that: 

J  J 


di\'<d(t¥< . <d& 


(10) 


On  Ek  to  compute  the  preceding  Xk  is  not  difficult  (see  the  following 
section),  and  moreover  this  kind  of  Xk  still  has  a  following  extremely 
important  characteristic . 

Theorem:  Under  stable  convergent  conditions,  if  Xk  is  the 
eigen  direction  of  system  (7)  on  Ek,  then: 


(11) 


In  which  q.  is  the  same  as  in  formula  (5) 
J  \ 

to  literature  £1j . 


(12) 

For  proof  please  refer 


III.  ALGORITHMS 


In  n  dimensional  space,  assuming  to  take  the  beginning  nxp  mat¬ 


rix  S  =(s. 

O  1 


(o) 


......  s^0^),  later  we  use  the  following  iteration 

form  to  progressively  produce  a  new  Sk: 
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1.  Evaluate  the  equations,  to  obtain  Zk: 

KZ.-S*. » 

2.  Gram-Schmidt  perpendicularization,  to  obtain  Yk  and  Rk: 

MZk-YkR* 

R>  ia 

on  the  triangle. 

3.  Form  given  effect  matrix  Bk: 

Bt~RkR$ 

4.  Evaluate  the  eigenproblem  of  given  effect  matrix  Bk: 

CD*  is  the  diagonal 


5. 

Form  tne  next  necessary 

step  Sk: 

S.-V.Q* 

(13) 

If 

in  the  above  algorithms 

we  make : 

S,-MX, 

(14) 

then  it  is  not  hard  to  prove  that  Xk  is  the  eigen  direction  in  the 
Ek  subspace  in  system  (7).  In  fact, 

X;3fX,-<tfYlAr‘Y*Q*-/, 

Therefore,  by  the  above  theorem  we  know  that  the  convergence  rate 
of  algorithm  (13)  is  (11)  and  (12).  In  fav.t,  formulas  (11)  and  (12) 
are  the  convergence  rate  of  the  preceding  simultaneous  iteration 
that  each  algorithm  shares. 

The  following  firstly  explains  some  problems  worthy  to  point  out 
of  algorithm  (13)  and  of  aspects  realized  in  the  computer.  Later  we 
use  a  flow  chart  (see  Fig.  1)  to  describe  the  complete  computation 
process . 

1 .  Firstly,  p  test  vestors  are  selected  for  use  through  a  sto- 
chastic  process  'ty  ,  or  selected  through  the  Bathe  method  ^  . 


(15) 

(16) 


2.  Iteration  formula  2  is  one  which  is  brought  about  by 
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the  Gram-Schmidt  perpendicular  process.  In  fact,  if  we  make: 


Zk—XkR*( Gram- Schmidt  perpendicular)  (17) 

In  which  X]£MXk=Ip,  Rk  is  the  upper  triangle  matrix.  The  two  sides 
of  the  above  formulas  simultaneously  the  remainder  takes  M  and  makes 
Yk=MXk,  thus  it  is  not  difficult  to  prove: 

YJM"Y»-I,  (18) 

The  computation  of  the  related  Yk  is  shown  in  Fig.  1 . 


3.  Because  algorithm  (13)  can  use  the  reduce  dimension  treatment 

technique,  therefore  while  some  certain  test  vectors  in  a  short  time 

satisfy  the  convergence  conditions,  we  will  not  again  let  them  add 

later  iteration.  But  due  to  the  reason  for  computation  error,  in  the 

later  iteration,  we  still  -must  on  the  other  test  vectors  "purify"  and 

shed  these  already  converged  vectors.  Additionally,  we  yet  must  point 

out  the  Bk  exponent  of  algorithm  (13)  will  continually  reduce  with 

the  increase  of  already  converged  numbers  of  test  vectors.  The  upper 

triangle  matrix  Rg  in  Fig.  1  is  the  left  lower  angle  p-q  exponent 

matrix  of  the  p  exponent  upper  triangle  matrix  R=(r..O.  Its  g  is 

ij 

an  already  converged  eigenface  number. 


4.  After  iteration  step  1,  when  the  obtained  q  eigenface  aj.  e 

( 1  )  *  • 

completely  converged,  the  diagonal  element  d^  in  D1  is  simply  the 

similar  eigen  value  ,  and  the  series  sP  ^  (j  «1 , 2, ......  ,g)  of  SI 

( 1  )  J 

still  is  not  the  similar  vector  x .  of  the  eigen  value  x.,  but  it 

J  J 

is  not  difficult  to  see  that  the  two  have  the  following  relationship: 


•r-Mxr  </- 1, . ,  o 


(19) 


In  order  to  obtain  xjl )  we  still  must  make  the  following  few  operations 

(20) 


(/  -u . . 


(21  ) 
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Generally,  test  vector  number  p  is  determined  by  the  following 
equation : 


P  “ min {9  +  8,  29} 


(22) 


5.  Criterion  for  convergence  can  be  refered  to  in  literature 
^2)and(8) .  The  criterion  for  the  proposed  effect  in  literature 


(23) 


In  which  c is  the  allowable  error  of  the  eigen  vector. 

6.  X,Y,Z,  and  S  in  the  flow  chart  can  apply  the  identical  storing 
up  position. 

7.  Because  various  algorithms  of  simultaineous  iteration  all 
can  apply  multiple  sum  formula  or  primary  point  displacement  to  con¬ 
duct  acceleration,  it  is  not  recounted  again  it  the  flow  chart.  The 
interested  reader  can  refer  to  material  1 ,  2,  and  6  . 

8.  Because  the  M  formation  in  the  iteration  process  remains  in¬ 
variant,  the  proposed  effect  eliminates  the  compressed  storage  method 
of  the  whole  zero  element. 

IV.  METHOD  COMPARISON 

In  this  section,  we  will  compare  the  algorithm  in  the  text  with 
several  other  frequently  used  algorithms.  Before  we  compare,  we  first 
will  explain  several  points: 

1 .  Assuming  bv  and  b  separately  to  be  the  half  tape  width  of 
j£  m 

K  and  M.  Although  in  the  algorithm  of  the  unresolved  M,  M  can  pro¬ 
gram  storage  according  to  compression,  but  in  order  +'  still  make  it 
convenient  we  still  use  bffl  to  count  operational  quantity. 


FIG.  Is  FLOW  DIAGRAM  FOR  THE 
PROPOSED  ALGORITH 
KEY:  (a)  Selected  first  nxp 
matrix;  (b)  K  formation  triangle 
solution;  (c)  Already  converged 
eigenface  number  g  =  0;  (d)  iter¬ 
ation  frequency;  (e)  When  j>1  ; 
(f)  (p-g  exponent  to  the  angle 
formation);  (g)  Convergence 
determination,  revised  g. 


2.  In  the  algorithm  in  which  we  can  use  the  "reduce  dimension" 
treatment,  when  counting  each  step  of  the  computational  quantity, 
we  use  the  equal  value  method.  If  we  assume  a  convergence  that  al¬ 
ready  has  i  eigenface,  at  this  moment  the  computational  quantity  for 
each  step  of  iteration  is  T(p-i),  then  the  average  computational  quan¬ 
tity  for  each  step  of  iteration  is: 
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- "  r-  » sv’-.w- 


Algoritha 

1 

Me CnW" 

Proposed 

Iteration 

Formula 

AZt-r,., 

*.-*?*. 

QfaQfDi' 

rt.Z,Q,Dk 

.  *••*«*» 
<z{z,-/*, 

*•-*.*1 

>W.0»  . 

KZ,.MX,., 

z,-zImz, 

qU,q,.d,' 

X,  ■  z,q,o, 

xz«.  mx,., 

z,-z,a, 

(  Zj.MZ,  .  /X, 

*••*»*{ 

Q»z,Q,  *  o,1 

x,  -  Z,Q, 

KZ,.X,., 

K,-MZ, 

KfXl.tZ, 

MfZlr, 

|  mz,  -  r,«, 

*»•*,*« 
qJ*,o,  -  o;‘ 
s,.r,Q, 

! 

Resolution 

1  of  M 

Resolved 

*4»L1X 

A-L-'KL-T 

r-M 

Resolved 

Af  -  LI?  , 

AmflKL‘T 
'  YmtJX 

Gnresolvej 

Unresolved 

Unresolved 

1  Unresolved 

1 

1 

j 

Projecting 
Matrix  For 

unable  to  us 

able 
to  use 

unable 
to'  use  ■ 

able 
to  use 

1  unable  to 

1  _ use  _ 

able  to 
use 

_  Standard 

IQ 

Standard' 

Standard 

Standard 

Generalize^!  Standard 

* 

Computaticl 
Quantity  d 
Each  Iterd 
ation 

■ 

<*<9,  •  »_J  *  «»* 

t*p<t,  * 

|  •/<*,»*.>  *»*» 

TABLE  1 :  COMPARISOH  OF  SEVERAL  ALGORITHMS 
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pY.t(p-'  > 

__UL» - - -  (24) 

P 

3.  We  use  J(p)  to  represent  the  computational  quantity  that 
uses  the  whole  eigenproblem  for  evaluating  symmetrical  matrix  p. 

We  use  2J(p)  to  similarly  represent  the  computational  quantity  of 
corresponding  generalized  eigenproblems . 

Examples  of  comparative  conditions  of  several  algorithms  are 
given  in  Table  1 . 

From  the  table  we  can  clearly  obtain  the  following  conclusions: 

The  algorithm  proposed  in  this  text 

1 .  makes  resolution  with  the  first  order  M  formation  less  than 
the  algorithms  of  Rutushauer  and  Reinsch.  Moreover,  when  M  is  in 
accordance  with  the  compression  program  storage,  the  computational 
quantity  of  each  step  of  iteration  in  reality  must  be  less  than  the 
computational  quantity  of  Reinsch. 

2.  in  each  iteration  step,  is  less  than  the  algorithms  of 
McCormick  and  Nicolai  in  making  the  product  of  the  first  order  M 
formation  and  the  test  vector.. 

3.  is  comparable  to  Bathe's  algorithm,  except  that  it  can 
carry  out  the  reduce  dimension  treatment,  and  doesn't  need  each  step 
of  iteration  to  evaluate  the  generalized  eigenproblem  of  p. 

V.  NUMERIC  VALUE  TESTS 

In  order  to  test  the  dependability  and  efficiency  of  the  algo¬ 
rithm  in  this  text,  we  applied  FORTRAN  language  on  a  SIEMEMS-7760 
computer  and  did  numeric  value  test  operations  of  many  areas.  Now 
we  will  list  examples  of  related  conclusions. 


1 .  Computed  Example  1 :  Applying  the  Reisch  method  and  the 
method  of  this  text,  we  computed  the  eigenproblems  of  certain  complete 
aircraft.  The  main  point  of  this  computed  example  is  the  convergence 
method  that  is  used  to  compare  the  two  methods  (See  Table  2)  and  the 
conclusion  precision  (See  Table  3).  The  test  vectors  that  this  com¬ 
puted  example  used  were  25,  and  obtained  17  eigenfaces.  The  17  eigen 
values  that  were  required  only  had  differences  measured  in  the  tenth 
place  (on  the  tens  place  after  the  decimal  point),  and  limited  by 
length  were  not  again  listed  in  the  table. 


Iteration 

f-rsn  u  on  r».v 

1  i 

i 

2 

D 

n 

5 

i  « 

7 

t 

6 

.  i 

10 

Converg- 

I  ¥ 

■  text  1 

0 

0 

° 

< 

>  4 

| 

M 

IS 

17 

ence 

numb«£- 

Reinjcb  ; 

0 

0 

0 

2 

u 

15 

16  | 

; 

16 

u  ! 

TABLE  2:  COMPARISON  OF  REINSCH'S  ALGORITHM  WITH  THE  PROPOSED  ONE 

2.  Computed  Example  2:  Exponent  is  500;  half  tape  width  of 
K  and  M  is  301.  This  computation  is  mainly  used  to  assess  the  com¬ 
puter  efficiency  of  each  algorithm.  It  uses  13  test  vectors  and 
computes  5  eigenface  (See  Table  4). 


Eigenvector  - 

1 

|  2 

1 

»  ! 

4 

5 

1 

0 

r 

0 

This  text^o-4) 

0.0004 

; 

0.0143  ; 

i 

0.0550 

0.0062 

0.0020 

0.0402 

0.0004 

0.0007 

ReioackOO'4) 

0.0100 

0.1221 

|  0.0652 

0.0040 

0.0160 

0.0251 

0.0210 

0.0447 

TABLE  3:  COMPARISON  OF  REINSCH'S  ALGORITHM  WITH  THE  PROPOSED  ONE 
IN  EIGENVECTOR  ACCURACY. 


Algorithm 

Ratishaacr 

Rtiaack 

McConakk 

Nicolai 

Batha 

Text 

Time-seconds 

1100 

1 

1669 

Ml 

070 

700 

TOO 

TABLE  4:  COMPUTER  RUN  TIME  OF  SEVERAL  ALGORITHMS 
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3.  Computed  Example  3:  Separately  uses  Bathe's  algorithm  and 
the  proposed  one  to  compute  three  eigenproblems .  This  example  was 
computed  in  order  to  confirm  that  this  text's  algorithm  has  higher 
efficiency  than  Bathe's  when  the  required  eigenface  are  comparatively 
many  (See  Table  6). 


Problem 

Exponent 

— HSTT - 

Tape  Width 

Test  Vector 
Number 

Required 
Eigenface  Numb 

1 

100 

s 

30 

25 

2 

200 

20 

IS 

25 

2 

10 

St 

30 

TABLE  5:  PARAMETERS  FOR  THREE  PROBLEMS 


Problem 

i 

2 

‘  3 

This  text-secpnds 

«s 

224 

551 

seconds 

194 

sts 

■ 

TABLE  6:  COMPARISON  OF  BATHE'S  ALGORITHM  WITH  THE  PROPOSED 
ONE  IN  COMPUTER  RUN  TIME 


VI.  CONCLUSION 

This  text  proposed  simultaneous  iteration  improved  algorithm, 
proved  its  convergency  rate,  discussed  the  problem  of  performing 
this  algorithm  in  the  computer,  and  made  comparisons  with  several 
frequently  used  simultaneous  iteration  algorithms.  We  can  conclude, 
the  scale  and  dimensions  of  system  (1)  and  the  tape  width  are  ex¬ 
ceedingly  large.  Comparison  of  the  proposed  algorithm's  computer 
efficiency  with  those  of  Rutishauer,  Reinsch,  McCormick,  and  Nicolai 
then  is  also  exceedingly  high,  and  the  required  eigenface  of  system  (1) 
are  many.  The  computer  efficiency  of  the  proposed  algorithm  compared 
with  Bathe's  algorithm  is  exceedingly  high. 

In  the  research  period  of  this  method,  we  obtained  Comrade  Guang 
De's  concrete  guidance  and  support,  and  obtained  the  help  of  Li  Guang- 
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quan.  In  this  we  deeply  express  our  gratitude. 
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SOME  PROBLEMS  OF  PLATE  BENDING  HYBRID  MODEL  WITH  SHEAR  EFFECT 


Wu  Changchun,  Chinese  University  of  Science  and  Technology 
Abstract : 

A  triangular  plate  bending  hybrid  element  is  constructed  according, t>  £ei»~ 
aer's  principle.  It  has  9  degrees  of  freedom,  and  takes  theJshear  iffntjat*  ■ 
account.  The  C*  formulation  of  the  thick  plate  is  extended  to  the  thin  plate 
and  the'shear  locking'phenomenon  can  be  avoided  through  imposing  the  discrete 
type  of  Kirchhoff's  constraints.lt  is  interesting  to  note  that  the  kind  of  dis¬ 
crete  constraints  can  be  directly  imposed  by  means  of  Hermite  interpolation*  , 
This  is  one  of  the  features  which  distinguish  the  paper  from  ref.  C  5  JandC  1  ). 
On  the  other  hand,  the  specific  matching  between  bending  and  shearing  is  discus-" 
sed  in  the  analysis  of  the  thick  plate.  It  is  shown  that  the  reasonable  m.hAi.t 
between  bending  and  shearing  complementary  energies  may  be  destroyed  by 
inhomogeneous  magnification  of  the  discretisation  errors.  Furthermore,  a  princi¬ 
ple  of  enerKy  regulation  is  pat  forward,  with  the  aid  of  which  this  matching 
problem  is  effectively  treated  so  that  the  unified  analysis  of  the  thick  and  this 
plates  is  achieved.  Finally,  numerical  examples  are  given. 


I.  FORWARD 

In  recent  years,  the  limited  primary  analysis  of  hard  plate 
with  shear  effect  and  the  general  primary  research  of  thick  plate 
construction  i3  receiving  more  and  more  serious  attention.  This  is 
no  doubt  because  the  deserved  question  for  discussion  has  extensive 
practical  value.  But  what  is  even  more  important  is  in  view  of  its 
theoretical  significance.  The  essence  of  this  kind  of  research  is 
the  problem  of  using  the  simplified  C#  sample  to  approach  the  C'  cat¬ 
egory. 

As  to  the  displacement  method,  the  key  of  constructing  commonly 
used  elements  lies  in  wanting  to  make  the  thick  plate  formula  indep¬ 
endently  selected  deflection  field  ®  and  corner  9i  able  to  automatic¬ 
ally  satisfy  the  zero  shear  strain  thin  plate  time  conition: 
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(1) 

Otherwise  the  so-called  shear  locking  phenomenon  will  be  able  to 
emerge  and  bring  about  an  ordinary  solution.  Various  limited  model 
problems  all  can  emerge  from  this  kind  of  condition  ^  . 

In  order  to  bring  about  the  thin  plate  constraint  (1),  we  can 
use  Yj  ,  and  their  derivatives  as  basic  variables  ^  ,  define  the 
independent  deflection  field  and  shear  strain  field,  and  the  corner 
thus  is  defined  by  formula  (1).  This  way,  when  plate  thickness  h-^0, 
shear  rigidity  tends  toward  oo,  and  naturally  we  have  the  result  of 
y.-*0.  Literature  ^3  and  4>  thus  acquire  very  great  success.  Their 

J 

inadequacy  is  that  they  gain  numerous  and  jumbled  units.  Their  rigid 
matrices  separately  are  20  and  36. 

As  to  Kirchhoff's  constraint  (1),  if  it  can't  be  strictly  sat¬ 
isfied,  then  the  discrete  type  of  Kirchhoff's  constraints  can  be  im¬ 
plemented.  In  literature  {^5^  the  best  selection  proposes  the  concept 
of  discrete  Kirchhoff  assumptions,  causing  formula  (1)  to  only  be  im¬ 
plemented  in  the  unit's  section  point  place,  thus  we  take  the  sample 
of  the  C°  type  and  successfully  apply  it  to  the  thin  plate  bending 
problem  of  the  C'  model.  But  because  this  method  must  omit  shearing 
deformation,  it  is  not  suitable  to  thick  plates. 

Rationally  utilizing  the  reduction  integral  method  is  helpful 
to  resolve  the  thick-thin  plate  problem.  Literature  (_6 }  proposes 
that  only  in  the  energy  function  of  matching  plates  can  the  shearing 
deformation  part,  if  to  the  integral  sum: 

.  •  J  J  (2) 

reducing  the  numeric  value  integral  exponent,  this  therefore  is  the 
"selective  reduction  integral  method".  In  this  way  we  retard  the 
rigidity  problem  of  the  unit  when  the  plate  is  thick,  and  the  integral 
sample  point  of  the  unit  implements  the  discrete  type  of  Kirchhoff 
constraint.  If  we  regard  the  shear  rigidity  a  in  formula  (2)  as  a 
penalty  factor  treatment,  then  the  selective  reduction  integral  method 
is  simply  the  penalty  function  method^.  To  make  this  kind  of 
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method  suitable  to  plates  of  different  thick-thin  proportions,  we 
must  select  a  different  numeric  value  integral  exponent  and  a  differ¬ 
ent  penalty  factor  value.  Clearly  this  is  an  experiential,  very 
strongly  improvised  process. 

This  text  attempts  to  avoid  the  displacement  method,  and  start¬ 
ing  out  from  the  thick  plate  mixed  variable  formula  construct  a  cur¬ 
rent  stress  hybrid  component  of  a  shearing  triangle.  In  this  area, 

(7-9) 

Cookv  has  done  forerunning  work,  systematically  analyzed  and  com¬ 
pared  this  type  of  numeric  value  performance  of  elements.  This  text 
purports  to  further  advance  the  work  in  this  area. 

II.  GENERAL  HYBRID  ELEMENT  OF  9  DEGREES  OF  FREEDOM-UH 

Here  we  will  directly  give  guidance  to  the  unit  rigidity  matrix 
of  the  element,  and  various  discussions  will  be  carried  out  in  each 
section  remaining  after.  The  Reissner  functional  equation  of  the 
shear  effect  is: 

«•-]*]’  CA##/8h+O/(0|+®i)  —Bi(Mii)  —  Bt(Qi) —pat'idA 

-  J  c*  C(®  —«)<?.+  (8„-9.)Af.+  (0/-  b'iMJida 

-  j*  c,  (  o.o> + M.0. + him,  0,)  da  ^  ^ ) 

Where  and  B^  separately  express  the  bending  complimentary  energy 
density  and  the  shearing  complimentary  energy  density  of  the  plate 
element.  For  the  significance  of  other  symbols  refer  to  literature  (J2^, 

Utilizing  the  balanced  equation  ^-Q^=0  eliminates  shearing 

stress,  and  makes  partial  integrals,  when  the  geometric  boundary  con¬ 
ditions  are  satisfied,  the  entire  functional  equation  of  the  dispersed 
system  can  be  expressed  as  : 

n**yi  {“  J^j* C(Af,w,+p)a>  +  BtiM „) + Bt(M,h,yidA 

+$9*  <3-*+  c.  ( o.m + M.0.+ fiw 0,) </* j 


(4) 


where  *>  is  the  weight  of  the  vector  of  the  element  bound¬ 

ary  point  outside  the  unit.  If  we  ensure  the  continuity  of  the 
element  common  boundary  point and-0r  ,  then  by  the  fixed  common  con¬ 
ditions  0-  we  can  obtain  the  crosswise  balanced  equation  of  the 
plate  element,  geometric  equations,  and  surface  mechanics  boundary 
conditions.  Aside  from  this,  we  still  will  obtain  the  balanced  condi- 
tion  of  the  element  common  boundary  pressure  point ^  '  .  Reflecting  this, 
we  can  apply  the  dividing  continuous  internal  force  field  to  formula 
(4). 


According  to  what  Fig.  1  expresses,  each  variable  referred  to  can 
define  the  following  test  solution  fields: 


The  linear  moment  of  force  field  in  the  element* 


M,  Af„DT-C4>MP> 

<0>-CP,  Pi  -  0»3r 


C  40- 


-  1 
0 

-  0 


x  y  o  o  o  o 
o  o  .  i  *  y  o 
oooooi 


o  0  ' 
o  o 
*  y  - 


(5) 


FIG.  1:  PARAMETERS  AND 
NODAL  VARIABLES  OF  UH 
MODEL 


•When  Mtj  obtains  linear  dispersion,  Mij,ij  in  formula  (4)  then  disolves. 
At  this  moment  formula  (4)  is  equivalent  to  the  improved  complimentary 
energies  function  formula  that  the  following  study  earliest  applies. 
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The  cubic  deflection  deflection  field  on  the  element  boundaries: 

{<»>  — C3,«  5ts  c5ai3T=— CSfD  <  ft  >  ^ 

and  {  & }  *■  C®i  9.i  9,i  Q**  9»*  8,*  oo,  0.,  0jr*J T 

Here  S„  expresses  the  a>  dispersment  of  the  element  i-j  side,  and 
it  is  taken  as: 

Sll—Hiflit+H ;/®/  +  H//(9/)/  +  H/»(8/)f  ^  rj  ^ 

Here  Hi j ,  Hij ,  etc.  are  the  cubic  Hermite  interpolation  value 
function.  The  linear  directional  rotational  angle  and  the  linear  torse 
angle  field  of  the  element  boundaries  separate  into: 

<8.>  -  C0.„  9.,,  9 CN J  { &  >  (  8 ) 

{8,>  *  (9/i*  9/**  9/»i)T*“  C^/3  ( & )  ( 9 ) 

By  formulas  (5)^(9)  we  yet  can  deduct  the  following  matrix 
formulas : 

The  directional  bending  matrix  and  torse  matrix  on  the  element 
boundaries  separately  are: 

*1  C<t»/D  <  P  >  (1Q) 

The  bending  complimentary  energy: 


BtW,,)  — 0  >rC<t>r  -gfrC CK* ) < 


1  —  y  0 

v  10 

0  0  2  (  1  +  *) 


If  we  remember  that  —  then  the  shearing  complimentary 

energy  expression  is: 

(13) 
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In  which  G=E/2(1+v)f  a  is  the  shearing  force  coefficient  (normally- 
taking  the  value  6/5).  The  constant  shearing  force  field  on  the 
element  boundaries: 


3...  SUF-CUCG'HPi 


CLl- 
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it 

mn 
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m„ 
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91 

m,, 
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lij  and  m. .  are  the  exterior  normal  line  direction  cosine  of  element 
i-j  side. 


With  (5)-^(14)  substitute  (4)  to  obtain: 


**' - |-{P>TC//»30}  +  {P}TC//».D{&>-{G}r{&>  (15) 

C//«D  -  J  J  <t>  ¥  C  C )  CO  )  + ■ 3  W  ^  ( 1 6 ) 

tHrf-f^aQ'yajHiJl  +  C&KNj  +  C&yCNjyds  (17) 

Jc;  <<o.>rc^+  «.>  wro+  (1g) 

Formula  (18)  is  the  equivalent  point  load  vector,  is  the  torse 

value  matrix  of  the  ©  field  inside  the  element.  This  internal  ®  field 
is  desirable  for  abic  dispersment  and  boundary  3  dispersment  to  be  mut¬ 
ually  alike,  but  according  to  the  principle  of  the  partial  sum  approach 
it  is  also  desirable  as  simplified  linear  dispersment. 

By  the  fixed  constant  conditions  of  formula  (15)  we  can  obtain  the 
formula  for  the  element  internal  force  and  the  balanced  formula: 


Eliminating  {p>  we  then  obtain  Cif){ &>  —  <(?}, in  which  the  element  rigid 
formation  (9x9): 


(20) 


c/o-c/^ycz/Mrc//..) 


III.  THE  DISCRETE  TYPE  OF  KIRCHHOFF  CONSTRAINT  AND  GENERALIZATION 
OF  MODEL  FORMULAS 

In  the  above  section,  the  <o  field  of  the  element  boundaries 
adopted  the  cubic  Hermite  interpolation  value  (7),  this  signifies 
in  point  i  we  have: 


The  detailed  dividing  yet  strengthening  with  the  network  of  this  kind 
of  Kirchhoff  constraints  is  equivalent  to  the  numeric  value  techniques 
in  literature  |  1  and  1  j .  Both  use  the  discrete  formula  to  implement 
thin  plate  constraints  (1),  only  however  here  the  treatment  is  even 
more  simplified  and  more  direct. 

Under  the  conditions  of  extremely  thin  plates,  the  UH  element 
also  cannot  emerge  to  solve  the  difficulty  of  the  locking  phenomenon 
or  the  numeric  value  abnormal  state.  A  living  example  of  numeric 
value  gives  the  thick  span  ratio  as  the  UH  solution  of  extremely  thin 
plate  of  1/10000.  We  can  see  that  the  numeric  value  is  stable,  the 
precision  is  good,  and  the  reduced  integral  method  under  extreme  con- 


-«*<=«*#>  4 

j  i 


FIG.  2:  DISTRIBUTION  OF 
DISPLACEMENT  SAMPLE  ON 
THE  ELEMENT  BOUNDARIES 


ditions  of  thin  plate  of  h^0  is  difficult  to  obtain  dependable  numer¬ 
ic  value  results. 
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We  want  to  analyze  the  use  of  thick  plate  with  the  stress  hybrid 
model.  We  must  satisfy  the  following  few  conditions: 

1.  The  computation  of  the  shearing  complimentary  energy; 

2.  The  independence  of  test  solution  9,  and  8  ; 

3.  The  matching  relationship  of  rational  winding-shearing. 

As  to  the  UH  element  of  this  text,  the  preceding  condition  "I" 
clearly  can  be  satisfied.  As  to  condition  "2",  by  Fig.  2  we  can  see, 
because  along  the  element  boundary  ab,  S '  takes  cubic  interpolation, 

5,  then  is  quadratic  distribution  and  test  solution  9,  thus  is  linear 
distribution.  In  this  way,  fully  at  point  a  and  b  constraint  (1)  is 
brought  on.  But  on  the  whole  element  boundary  section,  y,— 9,+5„=jM  , 
this  then  ensures  the  independence  of  the  UH  element  test  solution 
and  8.  The  hybrid  element  HTC  in  literature  (9)  in  the  same  way 
uses  the  cubic  deflection  field  (7);  but  the  above  element  is  not 
applicable  to  thick  plate  distribution  because  its  corresponding  dir¬ 
ection  rational  angle  definition  is  0,— S„  on  the  whole  element  boundary, 
therefore  condition  "2"  is  violated.  In  the  next  section  we  will 
carry  out  discussion  of  condition  "3". 

IV.  PRINCIPLES  OF  MATCHING  BETWEEN  BENDING  AND  SHEARING  AND  ADJUSTMENT 
OF  ENERGY 

The  problems  of  the  matching  relationship  of  the  hybrid  model 
exists  in  two  areas.  Firstly,  if'  a  stress  hybrid  element  increases 
the  number  of  the  internal  force  coefficient  B,  that  then  will  increase 
the  rigidity  of  the  element,  and  the  exponential  order  of  the  test 
solution  that  raises  the  displacement  then  will  increase  the  degree 
of  flexibility  of  the  element.  Therefore  the  problem  exists  of  match¬ 
ing  an  internal  force  and  exponential  order  of  displacement.  It  is 
common  knowledge  that  this  matching  must  satisfy  the  so-called  order 
condition  otherwise  mobil  deformation  will  emerge  and  thus  zero 
energy  model  state. 

As  to  thick  plates,  the  problem  will  become  even  more  complicated. 
Equation  (16)  can  be  expressed  as  : 
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(22) 


T 


-gfraHJ+ahKHj) 


In  which 


CH.D 


J  Jc<t>)TCCDC4>]J>4  ft  CH*>— 


Separate  and  complimentary  energies  and  are  mutually  correspond¬ 
ing.  Remember  the  discrete  errors  of  these  two  discrete  matrices  rising 
are  and  R^*  thus  the  equation  for  the  total  complimentary  energies 
of  the  element  can  be  expressed  as: 


JJ  (B^BJdA - -|frC(<P>rC 

^  *"  (23) 

+aV({9)TtHtU$)  +  Rt» 

Clearly,  the  shearing  complimentary  energy  R~  of  the  plate  element 

2  ^ 

will  be  amplified  h  fold  (corresponding  to  R^  ) .  Moreover  the  plate, 
the  more  the  discrete  solution  of  the  thick  standard  as  to  the  dev¬ 
iation  of  the  true  solution,  the  larger.  In  short,  the  rational  match¬ 
ing  relationship  of  the  bending  complimentary • energy  and  the  shearing 
complimentary  energy  will  meet  with  violation  due  to  unequal  ampli¬ 
fication  of  the  discrete  error.  Cook'  '  had  a  test  solution  for  the 
same  kind  of  displacement , changed  to  using  a  difforent  internal  force 
field  and  element  form  to  construct  four  kinds  of  hybrid  elements, 
among  which  only  "H5"  is  commonly  used  with  the  thick  plates,  the 
reason  being  the  imbalance  of  the  bending- shearing  match. 


Generally  speaking,  the  rigid  characteristics  of  the  preceding 
two  kinds  of  variations  of  matching  relationships  to  the  hybrid 
element  both  have  obvious  effects.  Therefore,  if  we  can  adjust  the 
coefficient  to  improve  the  element  rigisity,  we  can  also  go  through 
matching  adjustments  of  B^  and  B£  to  improve  the  element  rigidity. 
Clearly,  the  former  is  due  to  the  limitations  of  receiving  the  order 
condition,  the  adjustable  scope  of  p  is  very  small.  Moreover,  the 
adjusted  effect  is  not  continuous.  Contrarily,  the  latter  under  the 
conditions  of  unchanged  displacement  and  internal  force  test  solutions 
can  in  a  very  large  scope  continually  improve  the  rigidity  or  flex- 
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ibility  of  the  element.  Based  on  this  kind  of  consideration,  we  site 
the  following  energy  adjustment  principles. 

2 

As  to  the  error  sums  of  R1  and  ah  R2  of  equation  (23),  people 
have  much  difficulty  directly  controlling  their  proportion.  But  we 
can  indirectly  go  through  the  adjustment  to  the  shearing  compliment¬ 
ary  energy  to  equalize  their  influence.  A  simple  and  effective  method 
is  to  change  the  factor  .a  value  in  formula  (22).  If  we  make  0ra<6/5 
then  due  to  the  unequal  amplification  of  R^  and  R2  causes  the  excess 
shearing  complimentary  energy  to  be  reduced  and  consequently  eliminated. 
The  chosen  value  of  a.  can  be  determined  by  a  simple  computation,  to 
the  UH  element  we  can  take  <*  =  0.6,  if  we  make  a. =0,  then  there  is  no 
shearing  effect;  thus  the  UH  conveniently  deteriorates  into  the  uni¬ 
versal  thin  plate  hybrid  element. 

V.  COMPUTED  EXAMPLES 

To  examine  the  general  use  of  UH,  we  carried  out  comprehensive 
computations  of  the  square  plate  that  possesses  various  different 
thickness  span  proportions.  By  Fig.  3  we  can  clearly  show  the  mech¬ 
anism  of  the  above-mentioned  energy  adjustment  principles.  In  the 
figure  the  H  condition  a  solution  curve  of  UH  (dotted  line)  with 
a=1.2*0.6  and  gradually  draws  near  to  the  thick  plate  theory  solution 
curve.  Henceforth,  when  £>0,  the  curve  again  gradually  opens  into 
a  straight  line  and  parallel  or  mutual  with  the  thin  plate  theory 
solution.  Clearly,  the  transformation  of  °  is  not  sensitive  to  the 
effect  of  thin  plate  solution.  The  so-called  energy  adjustment  can 
only  be  significant  by  raising  the  thick  plate  analysis  of  the  con¬ 
trolled  effect  to  the  shearing  effect. 

Figure  3  simultaneously  draws  the  boundary  9  field  changing  to 
be  the  solution  curve  of  linear  distribution  time  (dot-dash  line). 

We  can  see  that  this  simplified  linear  •  test  solution  can  be  adopted 
to  the  thick  plate,  because  it  possesses  C#  continuity,  moreover  it 
ensures  the  independence  of  i  and  9,  .  But  it  does  not  completely 
satisfy  Kirchhoff’s  constraints,  therefore  it  is  not  applicable  to 


FIG.  3:  CENTER  DEFLECTION  CURVES  OF 
SIMPLY  SUPPORTED  SQUARE  PLATE  UNDER 
UNIFORM  LOAD  q 

KEY:  (a)  1/4  plate  6x6  mesh  size;  (b) 
Thick  plate  theoretical  solution;  (c) 
Thin  plate  theoretical  solution;  (d) 
Solution  when©  is  linear  distribution. 


thin  plates,  and  consequently  the  left  end  of  the  corresponding  sol¬ 
ution  curve  of  Fig.  3  follows  the  reduction  of  h  to  rapidly  deviate 
the  true  solution  and  distribution.  This  converse  example  further 
explains  the  necessity  in  UH  of  implementing  the  discrete  type  of 
Kirchhoff  constraints. 


(a) su  «*#$**« 

(b) sc  «*«**«* 
I  Cq  1 .  cu  ■*«**«* 


KEY;  (a)  Simply  supported 
plate'  under  uniform  load; 
(b)Simply  supported  plate 
under  concentrated  load;  (c) 
Consolidated  supported  plate 
under  uniform  load;  (d)Con- 
solidated  supported  plate 
under  concentrated  load. 


'?-■  •  *  -V  , 


FIG.  4;  CONVERGENCE  CURVES  OF  SQUARE 
PLATE  UNDER  UNIFORM  LOAD  q  AND  CENTRAL 
CONCENTRATED  LOAD  P  RESPECTIVELY  (h/s- 


=0.01,  N  is  tha  mesh  size,  R  is  the 
ratio  of  finite  element  solution  (UH) 
to  theoretical  solution  for  -  > 
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Tables  1-»4  and  Fig.  4  and  5  further  express  the  various  numeric 
value  energies  of  the  UH  element.  We  can  see  that  regardless  of 
displacement  or  internal  force  UH  uniformly  gives  satisfactory  re¬ 
sults  . 
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h/ s 

Thin  plate  theo¬ 
retic  solution 

UH 

literature 

3 

0.0001 

0.01160 

0.01141 

0.01 

0.01160 

0.01142 

0.01170 

0.05 
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0.20 
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0.25 
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0.01734 
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TABLE  2:  DEFLECTION  COEFFICIENT  rp  of  SIMPLY  SUPPORTED  SQUARE 
PLATE  UNDER  CENTRAL  CONCENTRATED  LOAD  P.  rp=WmaxD/Ps* ,  taking 
a=0.6,  v=0.3»  D  is  the  bending  stiffness. 
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TABLE  3:  CONVERGENCE  OF  THE  CENTER  DEFLECTION  OF  SQUARE  PLATE  WITH 
MESH  REFINEMENT  (UH:a=0.6,  h/s=0.01,  v=0.3) 
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A  LOCAL  STRAIN  FATIGUE  ANALYSIS  METHOD  AND  COMPUTER  PROGRAM 


Wu  Yisheng,  Institute  of  Mechanicsf  Chinese  Academy  of  Sciences 
Abstract : 

A  local  strain  fatigue  analysis  method  is  presented  for  evaluation  of  fati¬ 
gue  life.  First,  the  three  parameter  elements  of  load  increment*  strain  increment 
and  stress  increment  are  constructed  by  using  the  cyclic  stress— strain  curve  of  ■ 
the  material  and  the  cyclic  load  notched  strain  curve  of  notched  specimens.  Then, 
the  local  stress-strain  analysis  of  notched  specimens  under  complex  load  it 
made  by  means  of  these  elements  and  “availability  coefficient  matrix* given  by 
R.  M.  Wetzel.  The  damage  of  each  cycle  is  determined  on  the  basis  of  local 
strain  amplitude  and  correction  for  the  effect  of  mean  stress.  Last*  the  damage 
can  be  cumulated  according  to  Miner's  linear  cumulative  damage  theory  and 
the  life  can  be  evaluated. 

The  features  of  this  method  .consist  in  adoption  of  three  parameter  ele¬ 
ments  of  load  increment,  stress  increment  and  strain  increment!  one-step  transfor¬ 
mation  from  loadtime  histories  to  local  stress  and  strain-time  histories,  than 
the  computation  procedure  is  simplified.  ^  t. 

The  example  of  the  Cumulative  Fatigue  Damage  Division  of  thnSAE  Fati¬ 
gue  Design  and  Evaluation  Committee  has  been  computed  with  our  program.  The 
crack  formation  lives  of  the  notched  specimens  for  two  materials  under  vehicle 
transmission  load  have  been  evaluated  and  the  results  are  quite  consistent  with 
W.  R.  Brose's  results  and  experiments.  This  demonstrates  that  the  presented  nie«^ 
thod  and  program  are  simple,  reliable  and  rapid.  This  program  can  work  more 
than  50  peaks  per  second  and  may  be  available  for  engineering  evaluation  of 
the  crack  formation  life.  v-  -  *<* 

s  *  / 

I.  FORWARD 

Structure  fatigue  damage  analysis  and  life  estimated  computations 
are  always  a  problem  that  many  departments  of  industry  pay  close  at¬ 
tention  to.  Traditional  life  estimated  computations  have  adopted 
nominal  stress  methods.  If  the  life  is  computed  according  to  the 
nominal  stress-time  histories  that  the  component  receives  and  the  S-N 
curve  of  the  component  is  obtained  through  fatigue  tests,  it  is  gener¬ 
ally  necessary  to  pay  a  very  large  price.  This  kind  of  method  for 
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life  estimated  computations  cannot  correctly  reflect  the  damage 
process  of  the  stress  cumulative  position;  therefore  it  is  not 
accurate  enough. 

The  even  more  accurate  life  estimated  computation  method  that 
is  presently  used  is  the  total  life  of  the  component  divided  into 
crack  formation  and  crack  formation  expanding  into  two  phases.  Crack 
formation  phases  are  indicated  by  the  appearance  of  macroscopic,  ev¬ 
ident  cracks.  As  to  the  length  of  the  engineered  cracks,  this  is 
at  present  still  not  clearly  defined;  based  on  experience  this  can 
be  determined  as  2mm.  The  life  of  the  phase  of  crack  expansion  can 
use  the  method  of  break-crack  mechanics  based  on  material  of  crack 
expansion  characteristics  to  be  determined;  crack  formation  life  can 
be  determined  based  on  local  stress  and  strain  analysis  and  strain 
fatigue  material  of  the  stress  concentrated  position.  This  text  only 
touches  upon  crack  formation  life  computations.  In  the  construction 
of  excelling  material  and  referring  to  secure  life  designs,  crack 
formation  phases  are  extremely  important. 

Local  strain  fatigue  analysis  method  has  been  developed  by  for¬ 
eigners  in  the  last  ten  some  years,  has  already  begun  to  correspond 
to  departments  of  the  aeronautics  industry,  and  has  already  broadened 
to  the  automotive  industry  and  other  industry  departments  ^  .  This 
text  introduces  this  kind  of  method  and  its  computer  programs. 

II.  METHOD  OF  LOCAL  STRAIN  FATIGUE  ANALYSIS 

The  first  step  of  this  method  of  analysis  (also  the  key  step) 
is  the  transformation  of  load-time  histories  received  by  the  com¬ 
ponent  to  local  stress-strain  time  histories  of  the  stress  accumulation 
position.  Later  again  based  on  the  size  of  local  stress-strain  to 
compute  damage,  the  life  can  be  evaluated.  The  major  advantage  of  this 
method  is:  we  only  need  to  grasp  the  fatigue  test  data  of  the  material, 
then  we  can  forecast  the  life  of  each  type  of  component  under  various 
load  history  effects.  Additionally,  we  can  compute  damage  according 
to  the  sequence  of  added  load.  In  this  way  then,  we  consider  the 


inside  with  the  mutual  effects  between  loads.  Therefore  it  is  com¬ 
paratively  precise.  Under  the  conditions  of  stress  concentration 
position  with  model  deformation  and  cyclic  added  load,  due  to  the 
main  linearity  of  the  material  in  itself,  the  material  hardens  or 
softens  under  cyclic  load  increments,  and  the  mutual  effects  between 
loads  are  equal  factors,  causing  local  stress-strain  analysis  to  be¬ 
come  very  complex.  The  principle  says  that  using  the  model  flow 
theory,  the  method  that  uses  limited  elements  can  carry  out  cyclic 
load  increment  analysis.  But  at  this  moment  its  operation  weight  is 
too  large,  when  exceeding  the  cycle  a  few  times,  in  reality  there  is 
already  no  way  to  carry  it  out.  Therefore  we  must  look  for  a  suitable 
simple  method.  Based  on  large  amounts  of  experimental  research, 
present  local  strain  fatigue  analysis  uses  some  sample  assumptions 
that  follow; 

1.  Using  a  cyclic  stress-strain  curve  to  make  the  material's 
basic  stress-strain  relationship.  Because  the  hard  or  soft  phase 
of  the  cycle  of  the  material  in  the  fatigue  cycle  process  only  con¬ 
stitute  a  very  small  proportion;  after  going  through  a  determined 
cycle,  the  material  then  becomes  stable.  Therefore  we  can  use  a  cyclic 
stress-strain  curve  to  represent  the  interrelation  of  stress  size 
value  and  strain  size  value. 

2.  Cyclic  stress-strain  hysteresis  loop  and  cyclic  stress-strain 
curve  conditions  are  similar.  The  former  multiply  by  the  latter. 

3.  Load  cycle  and  local  stress-strain  cycles  correspond  one  by 
one.  If  the  load  completes  a  cycle,  the  local  stress-strain  also 
completes  a  cycle.  Therefore  we  can  construct  a  cyclic  load-strain 
curve,  and  this  curve  represents  the  relation  of  load  size  value  and 
local  strain  size  value. 

4.  The  material  possesses  memory.  If  after  the  stress-strain 
leaves  the  original  means  to  complete  a  cycle,  it  still  returns  to 
the  original  means. 


63 


(r.) 

T 

)  &«** 

a* 

tt*c 

KEY:  (a)  Load-time  history;  (b)  Load- 
notched  strain  transformation;  (c)  Cycle 
count;  (d)  Stress-strain  transformation; 
(e)  Accumulated  damage  computation;  (f) 
Fatigue  life  estimated  computations. 


(e)£ 


xfiattti-* 

r 


FIG.  1 :  THE  STEPS  OF 
LOCAL  STRAIN  FATIGUE 
ANALYSIS 


5.  Under  the  conditions  of  stress  concentration  position  in  local 
model  deformation,  because  the  model  deformation  area  is  surrounded 
by  the  surrounding  elastic  deformation  area,  therefore  is  similarly 
placed  in  strain  control  conditions  and  we  can  control  the  strain 
fatigue  experimental  data  to  carry  out  damage  analysis.  Regarding 
the  basic  principles  of  stress-strain  analysis,  please  refer  to  lit¬ 
erature  £2). 

General  steps  of  local  strain  analysis  are  shown  in  Fig.  1  ^ 
Firstly,  utilizing  cycle  load-strain  curve  with  the  load-time  histories 
received  by  the  component  to  transform  into  notched  strain-time  his¬ 
tories,  and  again,  utilizing  the  cyclic  stress-strain  curve  with 
notched  strain-time  histories  to  transform  into  notched  stress-time 
histories?  In  the  transformation  process  we  use  the  rainflow  count¬ 
ing  method  to  differentiate  cycles,  and  we  use  the  strain-life  curve 
to  compute  damage.  Later  according  to  Miner’s  linear  cumulative  dam¬ 
age  theory  to  cumulate  damage,  the  life  can  be  estimated. 


In  this  text  we  integrate  considerations  with  the  cyclic  load- 
strain  curve  and  the  cyclic  stress-strain  curve,  constructing  the 


•When  the  notch  is  in  single  axis  stress  condition  (For  example  a  thin 
plate  notch  specimen)  it  can  directly  cause  a  smooth  specimen  determined 
cyclic  stress-strain  curve;  when  the  notch  is  in  tri-axis  stress  cond¬ 
ition  (For  example  thick  plate  notch  specimen),  it  can  cause  limited 
element  computed  cyclic  notch  stress-strain  curve. 
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three  dimensional  cyclic  load-stress-strain  curve.  In  the  later 
computation  sequence,  we  use  the  three  parameter  element  of  load- 
strain-  stress  increments  one  step  transformation  from  load-time 
histories  to  local  stress- strain  time  histories.  From  this  step  of 
simplified  local  stress-strain  analysis,  it  shortens  the  time  of 
analysis  and  computation. 


FIG.  2:  LOCAL  STRESS-STRAIN  ANALYSIS  FOR 
CYCLIC  LOADING 


Now  using  Fig.  2..as  an  example  to  show  the  process  of  local 
stress-strain  analysis,  and  at  the  same  time  show  how  to  utilize 
the  above-made  assumptions.  Establishing  Fig.  2(b)  to  show  notch 
specimen  receiving  Fig.  2(a)  to  show  the  effect  of  load-time  historie 
In  order  to  conduct  local  stress-strain  analysis,  firstly  we  mus+-  go 
through  experiments  to  determine  cyclic  stress-strain  curve  of  the 
material,  as  the  dotted  line  in  Fig.  2(d)  shows,  and  go  through  ex¬ 
perimental  or  limited  element  computation  to  determine  the  cyclic 
load-strain  curve,  as  the  dotted  line  in  Fig.  2(c)  shows.  The  cor¬ 
responding  hysteresis  loop  uses  multiplie  principles  to  be  determined 
Analysis  from  the  load  peak  value  point  1  begins,  and  using  cyclic 
load-strain  curve  determination  and  the  largest  peak  value  load  cor- 
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responding  notched  strain,  if  point  1  in  Fig.  2(c)  assumes  the 
specimen  is  very  thin,  the  notch  similarly  is  in  single  axis  stress 
condition,  therefore  we  can  directly  use  the  cyclic  stress-strain 
curve  determination  and  the  notched  stress  corresponding  notched 
strain,  thus  point  1  in  Fig.  2(d).  Loading  from  point  1  removing  to 
point  2,  the  notch  strain  and  stress  separately  follow  the  hysteresis 
loop  transformation  that  sets  out  from  point  1  to  point  2  (the  point 
2  of  Fig.  2(c)  and  2(d)).  Loading  from  point  2  adding  to  point  3, 
notch  strain  and  stress  follow  the  hysteresis  loop  that  sets  out  from 
point  2  transforming  to  point  3  (the  point  3  of  Fig.  2(c)  and  2(d)). 
Loading  from  point  3  removing  to  2',  at  this  moment  loading  completes 
a  cycle,  notch  strain  and  stress  follow  the  hysteresis  loop  that  sets 
out  from  point  3  to  point  2'  (the  2'  of  Fig.  2(c)  and  2(d)).  At  this 
moment  stress  and  strain  also  complete  a  cycle.  When  continually  re¬ 
moving  loads  from  point  2*,  due  to  the  memory  of  the  material,  strain 
and  stress  break  and  return  to  the  hysteresis  loop  transformation  that 
sets  out  from  point  1  ,  goes  directly  to  point  4.,  and  the  later  process 
analogizes.  In  the  figure  using  2'  and  5'  to  express  the  points  both 
are  points  where  loading  and  stre  s-strain  complete  a  cycle.  In  local 
stress-strain  analysis  of  cyclic  ..oading,  correctly  determining  the 
completion  of  stress-3train  cycle,  and  correctly  selecting  the  memory 
of  stress-strain  means  to  reflect  the  materia] ,  both  are  extremely 
important . 

III.  WORKING  OUT  AND  COMPUTING  EXAMPLES  OF  COMPUTATION  TECHNIQUE 
AND  COMPUTATION  SEQUENCE 

The  analysis  of  local  stress  and  strain  is  very  complex.  The 
transformation  from  load-time  histories  to  local  stress  and  strain- 
time  histories  must  be  conducted  on  an  electronic  computer.  In  order 
to  evaluate  linear  problems  with  non-linear  problems  we  can  divide 
a  certain  amount  of  elements  with  the  cyclic  load-strain  curve  and 
the  cyclic  stress-strain  curve;  for  now  each  element  is  represented 
by  a  straight  line.  In  order  to  effectively  utilize  these  units  to 
compute  stress  and  strain  transformation,  and  correctly  reflect  the 
material's  memory  and  differentiate  the  completion  of  the  cycle,  the 
specified  useable  coefficient  that  composes  each  element,  we  adopt 
the  effective  coefficient  matrix  that  Wetzel  established  This 
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text  introduces  how  attached  load  increments  on  each  element  stress- 
strain  increment  establish  three  parameter  elements  of  load  strain- 
stress,  and  utilize  these  elements  with  load  transformation  to  com¬ 
pute  transformation  of  strain-stress.  The  following  simply  explains 
this  method  of  computation  and  the  steps  of  computation. 

1 .  Intake  source  data 

(1)  Load-time  histories.  Load-time  histories  use  the  sequence 
delivery  of  peak  value  and  valley  value.  In  order  to  obtain  completed 
closed  load  and  stress-strain  cycles,  the  best  piece  of  each  cycle 

is  from  the  biggest  load  beginning  to  the  biggest  load  end,  or  is 
from  the  smallest  load  beginning  to  the  smallest  load  end. 

(2)  Cyclic  load-strain  curve.  This  curve  generally  can  be  ex¬ 
pressed  by  the  following  functional  equation: 


(1) 


Here  C^  ,  C 2»  and  d  are  the  curve  drafting  together  the  constant; 
pa  and  e.  separately  are  cycle  size  value  and  strain  size  value.  With 
C.j  ,  C2,  and  d  three  constants  are  delivered. 

(3)  Cycle  stress-strain  curve.  This  curve  generally  can  be  ex¬ 
pressed  with  the  following  functional  equation: 


(2) 


Here  E  is  elastic  model  weight;  k'  is  cycle  intensity  constant;  Oh  is 
stress  size  value.  With  E,  k',  and  n'  the  three  constants  are  delivered. 


(4)  Strain-life  curve.  This  curve  can  generally  be  expressed 
by  the  following  functional  equation: 


)•+«', (2^,) 


(3) 


Where  Nf  is  the  violated  cycle  number;  au' ,  b,  and  c  separately 
are  called  fatigue  intensity  coefficient,  fatigue  intensity  index, 
and  fatigue  ductibility  index.  They  all  are  constants  related  to 
the  material  quality.  It  is  accumulated  by  these  four  parameters. 


r  (kg) 


FIG.  3:  CYCLIC  LOAD-STRAIN  CURVE,  CYCLIC  STRESS- 
STRAIN  CURVE  AND  RESPECTIVE  HYSTERESIS  LOOP  REP¬ 
RESENTED  BY  TEN  ELEMENTS 


2.  Establish  the  three  parameter  element  of  load-strain 
increment-stress  increment.  By  the  computations  of  formula  (l)  and 
the  related  notch  strain  of  the  highest  peak  value  load,  taking  these 
strain  increment  values  as  the  largest  strain  values,  by  the  cyclic 
load-strain  curve  and  cyclic  stress-strain  curve  are  divided  into 
strain  increments  equal  to  the  certain  amount  of  elements.  Based  on 
formula  (1)  and  formula  (2)  using  alternate  solutions  to  find  the 
load  and  stress  numeric  value  that  corresponds  with  each  point  strain, 
this  text  adopts  the  Newton  alternative  method,  using  the  given  load, 
strain,  and  stress  increments  that  each  element  represents.  The  hys¬ 
teresis  loop  also  respectively  divides  into  the  same  kind  of  numbered 
elements,  and  the  increment  multiple  that  each  element  represents  is 
obtained.  Figure  3  expresses  the  cyclic  load-strain  curve,  the  cyclic 
stress-strain  curve,  and  the  respective  hysteresis  loop  represented 
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by  ten  elements. 


3.  Establish  an  effective  coefficient  matrix,  compute  the 
transformation  of  stress  and  strain.  According  to  the  following 
rules  determine  the  available  coefficient  of  each  element: 

(1)  If  we  begin  from  the  largest  peak  value,  then  not  all  the 
element  available  coefficients  in  the  beginning  uniformly  take  +1, 
otherwise  they  take  -1 . 

(2)  The  available  coefficients  of  removing  load  time  elements 

by  +1  change  to  -1 ,  with  the  computed  load  descending.  When  removing 
load  we  only  can  use  the  elements  where  the  original  available  co¬ 
efficients  are  +1  ;  and  the  elements  where  the  original  available 
coefficients  are  -1  are  "unavailable  elements". 

(3)  The  available  coefficients  of  added  load  time  elements  by 
-1  change  to  +1,  with  the  computed  load  rising.  When  adding  load  we 
can  only  use  the  elements  where  the  original  available  coefficients 
are  -1;  the  elements  where  the  original  available  coefficients  are 
+1  are  "unavailable  elements". 

(4)  The  elements  are  each  used  once  and  the  symbols  change  once. 

(5)  Regardless  of  added  load  or  removed  load,  always  start  util¬ 
ization  from  the  first  element. 

In  the  computation  process,  whenever  the  condition  of  unavail¬ 
able  elements  is  met  with,  then  expressing  load  and  stress  and  strain 
completes  a  cycle.  The  emergence  of  unavailable  elements  becomes  the 
key  to  using  the  effective  coefficient  method  to  distinguish  stress- 
strain  cycles.  It  raises  the  identical  effect  with  the  rainflow  count¬ 
ing  method. 

4.  Compute  damage.  The  damage  that  each  stress-strain  cycle 
creates  can  be  computed  according  to  the  following  formulas: 
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1  / N 2  (ej/tj1"  e„>e„  ( 4 ) 


In  the  above  formulas  (Ta  is  stress  size  value;  (To  is  average  stress; 
while  t.,'  and  separately  are  elastic  strain  and  model  strain. 

They  can  be  computed  according  to  the  following  formulas: 

E  ((>) 

e,.-(a./A/),/*/  (7) 

5.  Accumulate  damage.  If  a  phase  of  load-time  histories  include 
n  stress-strain  cycles,  then  total  damage  can  be  accumulated  accord¬ 
ing  to  Miner's  linear  accumulative  damage  theory: 

II 

DAMSUM**  \/N,i  (8) 

1*1 

Here  1/Nfi  is  damage  that  the  i  stress-strain  cycle  has  created. 

If  this  phase  of  load-time  histories  is  repeatedly  exerted  on  the 
specimen,  then  the  violated  number  is: 

BTF-l/DAMSUM  (9) 

The  appendix  of  this  text  gives  the  entire  text  of  the  pro¬ 
cedure  that  is  drawn  up  using  Fortran  language  (omitted  here).  This 
procedure  is  the  improved  form  on  the  basis  of  literature  [_3)  •  What 
is  major  is  adopting  the  three  parameter  element  of  load-stress- 
strain  increments,  simplify  the  steps  of  computation,  and  increase  the 
automatic  dividing  procedure  of  the  element,  greatly  reducing  mach¬ 
ine  operation  time.  To  use  this  computation  procedure  refer  to  the 
examples  given  in  literature  ^5^.  Figure  2  shows  the  load-time  his¬ 
tories  that  the  notched  specimen  receives  constituted  of  1709  peak 
value  load.  The  largest  loads  separately  are  7265>3633»  and  1592  kg. 


70 


Specimens  were  manufactured  of  copper  by  the  trademarks  Man-Ten 
and  RQC-100.  Table  2  and  Fig.  4  and  5  we  can  3ee  that  the  results 
of  this  computation  and  the  results  of  computations  in  literature 
(3]  are  perfectly  close.  Moreover  they  coincido  with  experimental 
results.  Table  3  gives  computation  time.  Using  the  procedure  drawn 
up  in  the  text,  on  the  average  it  can  compute  50  peak  values  per 
second.  The  level  of  foreign  reports  is  computing  30  peak  values 
per  second^  . 
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TABLE  1 :  THE  RELEVANT-PARAMETERS  OF  SPECIMENS  AND  MATERIALS 
E,  k',  o£  unit:  kg/mm  ;  c^ ,  c2  unit:  kg. 


T-OT.loa-iv£omparison|  Violated 

load  .  '"CWfWsrife-  jCdmput 


Number 


suits  (thi-s- 


;er  re- 


eults  (3-) 


Experimental., 
rooulta  -(50 — 


72*5 


3033 


1IU 


text) 


15.2 


MM 


4900. 2 


154.1 


(009.0 


1.4 
12.  ( 
12.3 


420.9 

1(4.0 

74.0 


27SS.0 

4270.0 

1009.0 


TABLE  2a:  THE  COMPARISON  BETWEEN  COMPUTER  RESULTS  AND  TEST 
RESULTS  OF  MAN-TEN  NOTCHED  SPECIMENS 
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TABLE  2b:  THE  COMPARISON  BETWEEN  COMPUTER  RESULTS  AND  TEST 
RESULTS  OF  RQC-100  NOTCHED  SPECIMEN 
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TABLE  3:  TIMES  TAKEN  FOR  COMPUTING  1709  PEAKS 


BTF(violated  number) 


FIG.  4:  THE  COMPARISON  BETWEEN  COMPUTER 
RESULTS  AND  TEST  RESULTS  OF  MAN-TEN 
NOTCHED  SPECIMENS 

KEY:  (a)  kg;  (b)  Computer  results  (this  text); 
(c)  Computer  results  (3);  (d)  Experi¬ 
mental  results  (5). 


BTF(violated  number) 

FIG.  5:  THE  COMPARISON  BETWEEN  COMPUTER 
RESULTS  AND  TEST  RESULTS  OF  RQC-IOO 
NOTCHED  SPECIMENS 

KEY:  (a)  kg;  (b)  Computer  results  (this 
text);  (c)  Computer  results  (3);  (d) 
Experimental  results  (5). 


IV.  CONCLUDING  REMARKS 


1.  The  computer  results  make  clear,  the  life  that  the  computer 
obtains,  using  the  local  strain  fatigue  analysis  method  and  exper¬ 
imental  results  are  comparatively  close.  Using  this  method  to  com¬ 
pute  the  life  we  only  must  basically  test  the  data  of  the  material, 
and  therefore  it  is  more  advantageous  than  the  nominal  stress  method. 

2.  The  local  strain  fatigue  analysis  can  conveniently  compute 
accomplishments  on  the  machine.  The  procedure  drawn  up  in  this 
text  adopts  the  three  parameter  element  of  load  increment-strain 
increment-stress  increment.  Computation  is  rapid,  simple,  de¬ 
pendable,  and  can  be  used  to  estimate  computations  of  notched  spec¬ 
imen  crack  formation  life. 
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SHOCK  WAVE  BOUNDARY  LAYER  INTERACTION  IN  COMPRESSOR  CASCADES 


Yu  Shen,  Institute  of  Engineering  Thermophysics,  Chinese  Academy 
of  Science 


Abstract : 


Shock  *««  boundary  layer  interaction  in  compressor  cascades  is  an  ex¬ 
tremely  complicated  problem.  In  case  of  the  existence  of  the  separation,  it  is  on* 
of  the  most  important  factors  which  determine  the  performances  of  transonic 
compressor.  .  •  . 

However,  no  thoronch  investigation  has  been  made  on  shock  wave  boun¬ 
dary  layer  interaction  in  compressor  cascades  and  only  a  few  papers  on  this 
topic  have  been  published. 

It  is  shown  by' calculation  and  analysis  that  the  main  form  of  the  intera¬ 
ction  in  compressor  cascades  is  the  interaction  between  shock  wave  and  turbu¬ 
lent  boundary  layer  in  channels. 

The  analysis  has  made  clear  that  the  separation  criterion  proposed  by 
Pearcey  is  an  empirical  criterion  obtained  from  wind  tunnel  tests  of  a  RAE102 
airfoil.  The  conditions  are  quite  different  from  the  restrained  channel  flow  in 
compressor  cascades  with  strong  adverse  pressure  gradient.  As  for  compressor 
cascades,  the  boundary  layer  grows  rapidly  and  separates  at  lower  Mach  number 
upstream  of  the  shock  wave,  the  length  of  separated  region  and  the  location 
of  vortex  sheet  are  all  different  from  those  in  the  case  of  Pearcey' s  single  air¬ 
foil.  The  conclusion  has  been  drawn  that  it  is  inappropriate  to  apply  Pearcey' $ 
separation  criterion  directly  to  shock  wave  boundary  layer  interaction  in  com¬ 
pressor  cascades. 

It  is  proposed  that  the  separation  criterion  for  shock  wave  boundary  layer 
interaction  in  compressor  cascades  be  in  the  form  f  (Mi,  Pi/ Pv.t.t  Pr.t./ Pit  Re.) 
—  0,  which  involves  the  effects  of  pressure  gradients  in  front  of  and  behind 
the  shock  wave. 


I.  INTRODUCTION 


Shock  wave  boundary  layer  interaction  rising  separation  is 
one  of  the  major  factors  influencing  transonic  compressor  per¬ 
formance.  Shock  wave  boundary  layer  interaction  in  compressor 
cascades  is  an  extremely  complicated  problem,  however,  up  until 
now  there  still  has  been  no  deep  research  carried  out  on  it,  and 
very  few  papers  have  been  published.  Moreover,  it  is  simple,  that 
among  them  some  viewpoints  are  not  comprehensive.  Therefore,  we 
very  much  need  to  carry  out  systematic  and  deep  research  on  this 
very  important  topic. 

II.  THE  FORM  IN  CASCADE  SHOCK  WAVE  BOUNDARY  LAYER  INTERACTION 

Shock  wave  boundary  layer  interaction  has  many  kinds  of  forms. 
Straight  shockwaves  or  slanted  shockwaves,  laminary  boundary  lay¬ 
ers  or  turbulent  boundary  layers,  different  wall  surface  forms,  flow 
restraints  or  not,  shock  wave  strength  or  weakness,  boundary  layer 
separation  or  not?  Under  the  above  different  conditions,  the  forms 
of  shock  wave  boundary  layer  interaction  all  are  different.  Main¬ 
stream  intensity  is  different  conditions  of  high  supersonic  speed 
and  low  supersonic  speed.  The  same  is  low  supersonic  speed,  and 

the  forms  of  interaction  of  the  M  number  less  than  1.4  and  great- 

(l) 

er  than  1.4  can  also  be  different  '  .  Therefore  there  are  many  kinds 
of  forms  of  shock  wave  boundary  layer  interaction.  For  thirty  some 
years  research  work  on  shock  wave  boundary  layer  interaction  has  been 
carried  out  as  fully  as  possible,  but  it  still  appears  very  insuf¬ 
ficient  and  scattered. 

As  to  compression  cascades  we  can  say,  shock  wave  outcome  is 
either  with  laminary  boundary  layer  interaction  or  with  turbulent 
boundary  layer  interaction:  this  is  the  first  problem  we  must  re¬ 
solve.  The  two  have  great  differences  0. 

On  the  surface  of  the  cascade  blade,  the  boundary  layer  is  at 
the  outset  laminary  ;  after  going  through  a  distance  it  turns  into 
turbulent  boundary  layer.  The  length  of  this  laminary  area  can  be 
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long  or  short,  depending  on  Reynolds  number  Re  and  the  degree  of 
turbulence  of  the  flow  Tu.  But  from  the  front  blade  origin  all 
have  laminar  flow  area.  If  the  entrance  orifice  of  the  force  one 
compressor  (if  there  is  a  front  fan,  then  it  is  the  force  one  fan) 
the  degree  of  turbulence  generally  is  not  large,  and  possibly  as 
small  as  1$  or  0 . 5  ?  *•  ,  the  force  one  rear  motor  is  2?;  the  rear 

forces  are  5*^6 %  ,  even  reaching  approximately  10?.  This  all  is 

data  under  not  stalled  conditions.  Here  Tu  is  referring  to  V  7*/vr 
v  is  average  flow  velocity;  u'  is  velocity  fluctuation  value  of 
average  airflow  direction. 


Schlichting  '  ^  points  out  that  in  view  of  the  Reynolds  number 
Rec  the  blade  hypotenuse  is  approximately  10'’  to  10^;  under  high 
altitude  conditions  Re^  is  smaller.  The  author  chose  some  parameters 
to  test  compute  Re  .  The  results  that  were  obtained  are  listed  in 
Table  1.  They  are  approximately  3x10?  to  3x10°,  similar  to  Schlich- 
ting's  results. 


0) 


What  Kiock researched  was  the  position  of  the  transition 
point  under  conditions  of  low  speed  and  small  Reynolds  number.  When 
the  degree  of  turbulence  is  low,  from  the  blade  front  edge  rising 
directly  to  50?  hypotenuse  close  by  within  this  length  all  is  lamin- 
ary  boundary  layers,  as  seen  in  Fig.  1. 
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FIG.  1;  EFFECT  OF  Rec  AND  Tu 

ON  THE  TRANSITION  POINT 
KEY:  (a)  Laminar  flow;  (b) 
Separation;  (c)  Turbulent  flow 
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The  density  of  the  compressor  cascade  along  the  radial  changes, 
the  cascade  tip  density  of  different  compressors  also  is  different, 
generally  speaking  transonic  speed  compressor  cascade  tip  density 
is  constantly  between  1.05  and  1.35,  also  the  density  of  some  super¬ 
sonic  compressors  is  greater  than  the  above  scope.  As  to  the  com¬ 
monly  used  dual  circular  arc,  multiple  circular  arc,  advanced  com¬ 
pression,  and  other  forms  of  cascades,  based  on  the  selected  density, 
cascade  blade  installed  angle,  the  facing  angle,  and  the  M  number 
flow,  we  can  approximately  determine  the  position  of  the  shock  wave 
shooting  to  the  blade;  this  generally  is  the  back  half  of  the  blade. 
Therefore,  the  shock  wave  boundary  layer  interaction  in  the  cascade 
is  generally  the  interaction  of  the  shock  wave  and  the  turbulent 
boundary  layers. 

III.  SEPARATION  CRITERIA  OF  SHOCK  WAVE  BOUNDARY  LAYER  INTERACTION 

When  shock  wave  intensity  reaches  a  certain  degree,  it  can 
then  lead  to  boundary  layer  separation. 

On  flat  plates  when  the  M  number  of  straight  shock  waves  reach 
1.3,  it  then  can  lead  to  turbulent  boundary  layer  separation^*'. 

As  to  convex  surfaces,  bending  frequency  is  exceedingly  large,  sep¬ 
aration  area  is  exceedingly  long,  and  lower  pressure  is  exceedingly 
small . 

PearceyC^'  established  experiential  separation  criteria;  Boun¬ 
dary  layer  separation  occurs  on  the  wing  model  when  the  forward 
shock  wave  M  number  is  1.27  (Note:  this  is  not  the  free-flow  M 
number.  The  research  Pearcey  conducted  on  the  model  was  when  the 
free-flow  M  number  was  less  than  1,  airflow  on  the  wing  surface  ac¬ 
celerated  to  supersonic  speed.),  or  there  is  separation  when  the 
rising  static  pressure  of  the  shock  wave  rises  to  reach  Ps/p-j  is 
is  1.4.  Here  p^  is  the  static  pressure  of  the  forward  shock  wave; 
p  is  the  static  pressure  separation  point;  realistically  it  indicates 
the  static  pressure  when  there  is  separation  influence.  To  conduct 
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computations  according  to  the  static  pressure  when  there  is  sonic 
speed,  see  Fig.  2^'  . 


FIG.  2:  SURFACE  PRESSURE  DIS¬ 
TRIBUTION  BEFORE  AND  BEHIND 
THE  SEPARATION  POINT 


Up  till  now,  because  we  still  have  not  established  the  sep¬ 
aration  criteria  in  the  cascade  under  shock  wave  boundary  layer 
ineraction,  some  literature ^ ^  takes  the  Pearcey  separation  criteria 
and  directly  applies  it  on  the  cascade,  but  because  the  Pearcey 
criteria  are  the  experiential  criteria  that  the  tests  obtained  that 
were  conducted  on  the  element  wing  of  the  RAE  102  wing  model  at  the 
thickness  ratio  being  10^,  the  coming  flow  was  at  high  subsonic  speed, 
the  angle  of  attack  was  2® ,  acceleration  on  the  wing  surface  reach¬ 
ed  supersonic  speed,  and  passed  through  quasi-straight  shock  waves 
reducing  to  subsonic  speed.  Pearcey  conducted  tests  under  different 
forward  wave  M  numbers.  To  sum  it  up,  distribution  occurred  when 
p _/p1  was  1.4.  This  is  not  precise  enough.  In  reality  the  p  /p.. 
test  point  of  many  numbers  is  close  to  1.5.  Because  transonic  speed 
cascade  appropriately  moves  under  working  conditions  of  both  possible 
separation  and  possible  non-separation,  therefore  we  should  look  for 
more  precise  separation  criteria. 


80 


We  still  should  \ oint  out  that  Pearcey's  element  wing  test 
under  high  subsonic  speed  and  the  conditions  of  supersonic  speed 
cascade  are  very  different.  The  front  edge  of  the  supersonic  speed 
cascade  has  shock  waves,  the  entrance  orifice  of  various  cascade 
models  are  different,  the  degree  of  airflow  expansion  is  different, 
the  J  model  cascade  back  and  front  sections  are  even  and  straight, 
the  S  model  cascade  first  goes  through  a  phase  of  compression,  there¬ 
fore  the  development  of  the  groove  forward  shock  wave  boundary 
layer  and  Pearcey's  element  wing  condition  are  very  different. 

The  author  considers  that  we  especially  should  point  out  that 
in  the  cascade,  airflow  receives  constraints,  flows  inside  the 
groove,  the  counter  pressure  gradient  increases,  and  therefore  boun¬ 
dary  layer  increases  to  a  faster  velocity,  and  separation  can  then 
be  produced  under  smaller  forward  shock  wave  M  numbers;  when  there 
is  separation  the  shock  wave  forks  at  the  place  close  by  the  sonic 
speed  area.  The  length  and  vortex  of  the  separated  region  are 
different  from  Pearcey's  constrained  motion  of  the  element  wing. 

Therefore  the  author  considers  that  we  cannot  directly  apply 
Pearcey's  separation  criteria  to  compressor  aircraft  cascades. 

The  author  proposes,  as  to  compressor  cascade  flow,  we  should 
use  the  separation  criteria  of  the  included  interaction  forward 
shock  wave  pressure  gradient  (we  also  can  use  boundary  layer  thick¬ 
ness)  and  shock  wave  pressure  gradient.  Its  functional  relation¬ 
ship  should  be  f(Mr  PT.E./p2’  "ec)=0  (This  PL.E.’  PT.E.’ 

and  P2  separately  are  the  front  edge,  the  rear  edge,  and  the  stat- 
pressure  of  the  rear  shock  wave).  Because  the  forward  shock  wave 
M.j  is  the  primary  parameter  of  resolving  shock  wave  strength  and 
weakness,  it  naturally  should  be  one  of  the  parameters  of  separation 
criteria.  The  cascade  back  and  front  sections  have  many  kinds  of 
different  forms,  the  compression  level  to  the  airflow  is  different, 
and  this  has  a  fairly  large  influence  on  the  boundary  layer  increas¬ 
ing  velocity  and  separation.  Therefore  separation  criteria  should 
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include  (p^-p^  v  )  or  adopt  its  immeasurable  outline  form 
(p-j-p^  g  )/p-j  >  thus  1-(p^  j  /p-j);  also  we  can  use  the  interacting 
forward  boundary  layer  thickness  as  a  parameter.  The  tests  of 
Mateer  and  others'-  '  clearly  indicate  that  between  the  separation 
region  length  and  the  interacting  forward  boundary  layer  thickness 
there  is  a  direct  ratio  relationship,  but  calibrating  the  boundary 
layer  thickness  is  not  as  convenient  as  calibrating  pressure,  and 
therefore  we  still  use  p£  g  /P-]  as  a  parameter. 

The  major  effect  of  separation  to  pneumatic  performance  is 
causing  damage  to  increase.  From  the  practical  viewpoint  we  come 
to  see  that  what  is  serious  still  does  not  lie  in  completely  not 
producing  separation,  but  lies  in  avoiding  producing  comparatively 
large  separation  areas,  so  as  not  to  let  damage  markedly  increase. 
From  this  point  of  view  we  can  see  that  because  the  major  separation 
of  compressor  cascade  in  the  groove  is  turbulence  separation,  rear 
shock  wave  subsonic  expands  pressure,  counter  pressure  gradient 
is  comparatively  large,  and  has  a  very  large  influence  on  the  sep¬ 
aration  area.  Therefore  we  can  also  take  the  counter  pressure 
gradient  of  the  rear  shock  wave  as  one  of  the  parameters  of  separ¬ 
ation  criteria,  using  its  immeasurable  outline  form  (p^,  g 
According  to  the  test  results  of  Chapman  and  others  Til)  ,  as  to  tur¬ 
bulent  separation  we  can  say,  the  influence  of  the  Re  number  to  the 
beginning  separation  point  is  not  great,  but  the  Re  number  has  fairly 
great  influence  on  the  length  of  the  separation  area^  .  These  the 
foregoing  are  the  major  parameters  of  separation  criteria. 

Because  the  interaction  of  shock  wave  boundary  layers  is  per¬ 
fectly  complicated,  we  still  cannot  completely  go  through  theoretical 
analysis  to  obtain  concrete  separation  criteria.  But  conducting 
tests  in  a  supersonic  speed  cascade  wind  tunnel,  and  analysis,  we 
can  obtain  half  experiential,  comparatively  precise  separation  cri¬ 
teria. 
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IV.  CONCLUSION 


Shock  wave  boundary  layer  interaction  in  compressor  cascades 
is  an  extremely  complicated  problem.  Under  conditions  of  rising 
separation  it  is  one  of  the  major  factors  of  resolving  trans¬ 
sonic  pressure  performance. 

The  author  has  gone  through  simple  computations  and  analysis 
to  clearly  show  the  main  interaction  form  of  cascades  is  groove 
shock  wave  and  turbulent  boundary  layer  interaction. 

The  author  carried  out  analysis  of  Pearcey’s  separation  cri¬ 
teria.  Because  in  cascades,  boundary  layers  increase  to  a  faster 
velocity,  under  the  M  number  of  smaller  forward  shock  waves  we  thus 
can  have  separation.  Th  e  development  of  forward  shock  wave  bound¬ 
ary  layer,  separation  area  increase,  and  vortex  position  also  are 
different  from  Pearcey’s  element  wing  conditions.  Therefore,  the 
author  considers  that  it  is  not  suitable  to  directly  apply  Pearcey’s 
separation  criteria  to  compressor  cascades. 

Through  analysis,  the  author  proposed  to  use  separation  criteria 
on  shock  wave  boundary  layer  interaction  in  compressor  cascades  that 
includes  forward  shock  wave  pressure  gradient,  rear  shock  wave  pres¬ 
sure  gradient,  and  also  include  in  the  criteria  and  Re  numbers. 

The  form  is  f  (M1  ,  P-|/pL  •  PT  E  / P2*  Rec)=0.  We  can  conduct  tests 
and  analysis  in  a  supersonic  speed  cascade  wind  tunnel  to  obtain 
concrete  separation  criteria  for  shock  wave  boundary  layer  inter¬ 
action  in  compressor  cascades. 
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EXPERIMENTAL  STUDY  ON  DISCHARGE  AND  LOSS  COEFFICIENTS  OF  COM¬ 
BUSTOR  SWIRLERS 


Zhang  Chunxia,  Wuxi  Aircraft  Engine  Research  Institute 
Wang  Yulian,  Shenyang  Liming  Machinery  Company 

Abstract : 


The  blade-type  swirlers  have  been  used  in  many  air  engine,  combustors  to 
control  the  primary  airflow  of  flame  tube  and  establish  a  strong  recirculation 
zone.  The  losses  produced  by  the  swirler  depend  on  stagger  angle,  pitch/chord 
ratio  and  passage  area  of  blades. 

Some  test  results  of  the  swirlers  with  different  stagger  angles,  numbers  of 
blades  and  flow  areas  are  presented.  A  set  of  discharge  characteristics  and  rela¬ 
tions  between  loss  coefficient  and  swirl  number,  stagger  angle  or  blade  number 
are  provided.  In  addition,  the  correlation  of  loss  coefficient  to  discharge  coe¬ 
fficient  of  blade-type  swirlers  has  been  established.  It  is  found  out  that  the 
discharge  coefficient  varies  from  0.6  to  0.9  for  swirlers  with  different  stagger 
angles.  The  correlation^  ”31 .0(1 . 53s—  1  )  indicates  that  the  selection  of  swirl 
number  must  retain  the  loss  coefficient  at  an  acceptable  level  (aerodynamically). 

It  can  be  seen  from  the  test  results  that  the  blade  loss  is  evidently  smaller 
than  the  secondary  expansion  loss  for  a  given  combustor  swirler.  The  value  of 
the  blade  loss  is  about  18  percent  of  the  total  loss  for  a  straight  blade  swirler 
with  stagger  angle  varying  from  68'  to  ?0‘. 

The  tests  have  also  proved  that  the  twisted  blade  swirler  is  superior  to  the 
straight  one  from  the  aerodynamic  point  of  view.  But  its  manufacture  is  more 
complicated.  Therefore,  the  straight  blade  swirlers  (with  stagger  angle  45’ ) 
have  been  used  in  a  large  number  of  combustors. 


I.  TEST  PREPARATIONS 

Scheme  of  the  test  facility  is  shown  in  Fig.  1. 

Calibrated  air  discharge  of  the  test  facility  is  0.06-0.14  kg/seconds. 
Mach  number  of  swirlers  entrance  orifice  M=0.03~0.14. 

Inside  the  polymethyl  methacrylate  tube  the  diameter  of  the 
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FIG.  1:  SCHEME  OF  THE  TEST  FACILITY 
KEY:  (a)  Voltage  valve;  (b)  Air-liquid  separ¬ 

ator;  (c)  Voltage  valve;  (d)  Air  desiccator;  (e) 
Pressure  table;  (f)  Valve;  (g)  Rectifier;  (h)  Mer¬ 
cury  thermometer;  (i)  Opening  plate  forward  pres¬ 
sure  gauge;  (j)  Opening  plate;  (k)  Differential 
pressure  gaugge ;  (1)  Entrance  orifice  total 

pressure  gauge;  (m)  Entrance  orifice  moving  pres¬ 
sure  head  gauge;. (n)  Swirler  test  specimen; 

(o)  Polymethyl  methacrylate  tube;  (p)  Rectifier; 
(q,  r,  s)  Swirler  exhaust  orifice  static  pressure 
gauge;  (t)  Exhaust  orifice  total  pressure  gauge; 
Tu)  Exhaust  orifice  convergence  section. 


installed  swirlers  is  ♦  235mm,  taking  the  size  of  the  observance 
recirculation  zone  and  the  variations  of  its  edge  axial. 


In  front  of  the  opening . plate  flow  meter  and  behind  the  swirler 
(In  front  of  the  exit  orifice  pressure  calibration  section)  devices 
are  installed  that  have  many  openings  rectification  to  eliminate 
tone,  utilizing  the  U  form  tube  to  calibrate  each  pneumatic  parameter. 

To  eliminate  system  errors,  the  system  total  pressure  damage 
that  i.°  measured  when  there  are  no  swirlers  follows  the  variable 
value  of  the  entrance  orifice  airflow  parameters,  in  order  to  finally 
compute  when  arrang-ng  the  loss  coefficients  of  the  swirlers. 

When  experimenting,  air  temperature  is  within  2O-«30®  C. 
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II.  CONSTRUCTING  COEFFICIENTS  OF  TESTED  SWIRLERS 


The  main  experimental  specimen  is  the  straight  blade  swirler. 

For  the  purpose  of  comoarison,  experiments  were  also  conducted  on 
the  twisted  blade  swirler  of  the  combustion  chamber  of  certain  types 
of  engines.  Its  minute  details  of  construction  and  dimensions  are 
seen  in  Table  1  and  Fig.  2  and  3.  In  the  Table  each  dimension  uni¬ 
formly  is  a  realistic  measurement  value. 

To  conveniently  compare  the  swirl  intensity,  the  table  yet  lists 
uhe  swirl  number  S  of  each  swirler.  The  ratio  of  the  angular  mom¬ 
entum  and  the  axial  momentum  of  the  swirler  coefficient  exhaust  ori¬ 
fice  airflow  and  the  effective  diameter  product  of  the  swirler 
is  computed  according  to  the  following  formula: 


(D 


Where 

Rn  is  ring  radius  inside  the  swirler; 

R  is  ring  radius  outside  the  swirler; 
0  is  staggered  angle  of  the  blade. 


FIG.  2:  STRAIGHT  BLADE  SWIRLER 
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FIG.  3:  TWISTED  BLADE  SWIRLER 


TABLE  1  CONSTRUCTION  FACTORS  OF  TESTED  SWIRLERS 


'>Sv£arameter 

Number''^. 

Blade 

Staggered 

Angle 

B 

Blade 

Num¬ 

ber 

n 

Blade 

Thick 

ness 

mm 

b 

Swirler 
outside 
diameter 
mm  Ds 

Swirler 
inside 
diameter 
mm  d 

Swirler 

Exhaust 

Orifice 

Cross- 

section 

min  Fs 

Swirler 

number 

S 

A-1 

1  1 

5  ' 

1 .15 

40.0 

8.9374 

1.9698 

A- 2 

1  mi 

5 

1 .15 

m  ret*® 

40.0 

9.1964 

1.9108 

A  2-1 

6 

1.15 

Hr  ivl 

40.0 

8.0033 

2.1632 

A2-2 

6 

1.15 

mmt  irt '  Sij 

40.0 

8.6561 

1,9962 

A3-1 

69"  5' 

9 

1.15 

70.0 

40.0 

7.6956 

2.1068 

A3-  2 

68°  41  1 

9 

1.15 

70.0 

40.0 

7.8666 

2.0635 

B-1 

48*  36’ 

9 

1.15 

70.0 

40.25 

13.6985 

0.9133 

B-  2 

46" 56’ 

9 

1.20 

70.0 

40.25 

13.6233 

0.8615 

32-1 

52- 1  » 

10 

1 .20 

70.0 

40.30 

14-0542 

1  .0312 

B2-2 

49c  24' 

10 

1  .30 

70.0 

40.30 

14.8139 

0.9394 

B3-1 

52°  13' 

12 

1.20 

70.0 

39.90 

13.7514 

1.0387 

B3-2 

50*  21 ' 

12 

1.30 

70.0 

39.90 

14.5002 

0.9715 

C2-1 

57®  45 1 

6 

1.20 

70.0 

39.90 

12.7799 

1.2762 

C2-2 

57' 30' 

6 

1.20 

70.0 

40.25 

12.8886 

1.2639 

C3-1 

59*  9' 

a 

/ 

1 .20 

70.0 

40 . 30 

11 .7167 

1 .3431 

C3-2 

58" 26' 

9 

1  .20 

70.0 

39.30 

11.9757 

1 .3105 

81 1® 

79*  O 

5 

1 .20 

71 .0 

30.0 

4.30 

2.1800 

•  Twisted  blade  swirler 

0  Blade  staggered  angle  is  average  value. 


III.  TEST  METHODS 

1 .  Determination  of  Swirler  Discharge  Coefficient  ^s 

a.  According  to  calibrate  parameter  computations  opening  plate 
discharge  is 

b.  According  to  discharge  continuity,  G. =G  (G  is  air  discharge 

^  3  s  u 

going  through  swirler),  seeking  the  swirler  discharge  coefficient  rs. 
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2.  Determination  of  Loss  Coefficient  % 


According  to 


4l-AP?/(PJ-J>.) 


(2) 


Where 

4,  is  corresponding  to  the  airflow  loss  coefficient  of  the 
swirler  intake  orifice. 

AP?'  is  the  total  pressure  damage  going  through  the  swirler. 

Pi.— Pi  is  the  swirler  intake  orifice  moving  pressure  head. 

The  loss  coefficient  corresponding  to  the  swirler  exhaust  ori¬ 
fice  cross  section: 


4,-4, (PVP,)*  (3) 

Where 

Fs  is  swirler  exhaust  orifice  cross  section  area; 

F1  is  swirler  intake  orifice  flow  tube  cross  section  area. 

In  order  to  obtain  the  swirler  loss  coefficient  of  the  non¬ 
quadratic  expanded  damage  (thus  blade  damage  coefficient),  we  still 
adopt  the  following  test  programs: 

The  swirler  exhaust  orifice  directly  opens  wide  into  the  at¬ 
mosphere,  calibrates  the  sum  pressure  of  the  exhaust  orifice  passage¬ 
way  among  each  blade  of  the  swirler  following  all  along  the  circum¬ 
ference  and  all  along  the  diameter,  and  lastly  takes  its  arithmetic 
mean  to  compute  the  loss  coefficient. 

IV.  TEST  RESULTS  AND  DISCUSSION 

Test  results  are  shown  in  Fig.  4  through  Fig.  9. 

1.  Discharge  Coefficient  of  the  Swirler 

The  tests  clearly  show  that  the  primary  discharge  coefficient 
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of  the  swirler  is  related  to  the  constructed  model  that  was  adopted 
and  the  blade  stagger  angle. 

(1)  Comparing  the  test  results  of  twisted  blade  swirlers  and 
other  straight  blade  swirlers,  we  can  know  that  the  average  blade 
stagger  angle  of  the  one  most  fully  forward  reaches  79°  *  but  many 
straight  blade  swirlers  where  its  discharge  coefficient  is  even 
smaller  than  the  stagger  angle  are  all  large.  Estimating  this  is  due 
to  using  back  twisted  blade,  its  upper  area  side  edge  is  parallel 

to  the  airflow;  the  airflow  angle  of  attack  is  zero;  thus  we  obtain 
the  low  of  the  straight  blade  swirler  where  the  damage  is  greater  than 
the  angle  of  attack. 

(2)  The  same  constructed  model  swirler  discharge  main  coefficient 
is  related  to  the  stagger  angle  of  the  blade.  Under  the  conditions 

of  this  test,  the  approximate  range  is; 

a.  as  to  45d<9<r50‘,  ^s  =  0. 800. 90 ; 

b.  as  to  50*<  0  <60e  ,  ^  =  0.65^0.75; 

c.  as  to  60°'  P^703  ,^s  =  0.60~"0.70. 

As  to  the  twisted  blade  swirler  of  =79°  ,  ^s  =  0. 74*~0. 84. 

(3)  The  discharge  coefficient  follows  the  increase  and  reduction 
through  the  swirler  airflow  pressure  decline.  Within  the  range  of  the 
pressure  decline  of  the  tests  in  this  text  assume  a  linear  relationship. 
But  the  discharge  coefficient  reduces  very  slowly,  in  similar  com¬ 
putations  we  can  assume  each  kind  of  swirler  discharge  coefficient  is 

a  constant. 

2.  Characteristics  of  Swirler  Discharge 

As  in  Fig.  5,  each  characteristic  line  corresponds  to  a  p  and  n 
value.  Refer  to  the  figure. 

Model  coordinates  G—GyT\/F,P\ 


Where  Gs  is  air  discharge  of  each  second  through  the  swirler 
(kg/seconds ) ; 

T?  is  the  sum  temerature  of  the  swirler  intake  orifice  (K); 
l  2 

P?  is  the  sum  pressure  of  the  swirler  intake  orifice  (kg/cm  ); 

I  2 

Fs  is  the  circulation  area  of  the  swirler  (cm  ). 

Longitudinal  coordinates  are  p  —  (PJ— Pt)/P[ 

Where  p^  is  the  static  pressure  of  the  swirler  exhaust  orifice. 

The  checking  computations  of  the  flow  circuit  parameters  in  the 
combustion  chamber  can  apply  to  this  formula  curve,  by  5  check  p, 
to  obtain  P2q* 

3.  Loss  Coefficient  of  Swirler 

(1)  The  damage  of  the  swirler  is  composed  of  the  two  parts 

which  are  the  quadratic  expanded  damage  and  the  blade  damage.  Accord¬ 
ing  to  the  test  results  in  literature  ,  when  p  value  is  within 
65®~85° ,  the  blade  damage  of  the  straight  blade  swirler  approximately 
holds  23 %  of  the  total  damage.  Making  (  ,i,»  is  the  blade  damage 

coefficient,  4,  is  the  total  damage  coefficient),  the  tests  made  clear: 

a.  Blade  stagger  angle  is  bigger,  ^  is  bigger,  moreover  we  can 

take  60  as  the  dividing  1‘  ,  after  p<60° ,  *  unobviously  increases; 

b. Tha  *  value  distance  of  the  twisted  blade  is  smaller  than  that 
of  the  straight  blade; 

c.  As  to  the  straight  blade  swirler. 

When  9  =68*'-  70*  ,  i unif or„  =  0 . 1  8 

When  9  -50t~60”,  i  unifoPB  =  0 . 08. 

Because  the  total  pressure  calibration  error  of  the  blade  passage 
exhaust  orifice  is  fairly  large,  and  moreover  the  total  pressure  dis¬ 
tribution  in  the  straight  blade  exhaust  orifice  is  not  as  uniform  as 
that  of  the  twisted  blade  exhaust  orifice,  therefore  test  lesults 
only  are  approximate  values. 

(2)  The  relationship  between  the  swirler  loss  coefficient  and  the 
blade  stagger  angle,  blade  number,  and  swirl  number. 
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a.  With  blade  stagger  angle  increase  and  blade  number  increase, 
the  swirler  loss  coefficient  increases: 

*,-16.5(O.75mc‘0-  1  )  (0-45*  ~70*) 

b.  Swirl  number  S  is  an  important  parameter  when  the  swirler 
is  designed.  For  the  combustion  chamber  of  an  aircraft  motor  it  is 
necessary  to  make  3^0.6,  then  we  can  establish  a  stable  combustory 
recirculation  region.  But  with  the  increase  of  S,  the  swirler  ef¬ 
ficiency  drastically  reduces,  therefore  the  two  must  be  given  con¬ 
sideration.  Based  on  the  test  results,  this  text  establishes  the 
formula  for  the  relationship  between  the  swirler  number  and  the  loss 
coefficient : 


*,-3.10(1.535-1)  (5) 

4.  The  Relationship  Between  Swirler  Discharge  Coefficient  and 
Loss  Coefficient 

Through  linking  each  test  swirler  discharge  coefficient  and 
loss  coefficient,  we  obtain  the  curve  that  Fig.  9  expresses,  and 
can  express  it  as: 

1 

**'"  1+0.07^4,  (6) 

As  it  is  convenient  to  compare,  according  to  the  computation  of  sim- 

ilas  hydraulic  formula  |i  — - -  -  the  curve  is  shown  in  the  figure. 

1  +v  * 

V.  SUMMARY 

The  preliminary  results  that  this  test  obtained  are: 

1 )  Swirler  discharge  coefficients  are  different  than  blade  type 
and  stagger  angle,  regarding  the  straight  blade. 

(1)  H,-0.8— 0.9  (45*<0<5#*) 

(2)  0.86— 0.75  (SO*<0<«O*) 
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FIG.  4:  GENERAL  DISCHARGE 
CHARACTERISTICS  OF  SW-IRLERS 


FIG.  5:  VARIATION  OF  LOSS 
COEFFICIENT  WITH  STAGGER 
ANGLE 


(3)  n,*#.  60— 0.70  (60* <  0  <70* ) 

2)  Obtaining  a  formula  that  correlates  to  general  discharge 
characteristics  of  swirlers  of  different  stagger  angles  and  blade 
numbers . 

3)  The  proportion  of  swirler  blade  damage  to  total  damage  is 
very  small.  Under  conditions  of  this  test,  when  fl  =68^-70°  , 
♦unifor."0-18’  uhen  »*5Cf~60',  » uniform=0.08. 

4)  Between  swirler  loss  coefficients  and  blade  stagger  angle 
and  swirl  number  the  following  linking  formulas  are  established: 

♦, -10.5(0. 75sec*0-  1) 

*«-»!.  0(1.555-  1  ) 

5)  After  9^60°  ,  loss  coefficient  markedly  increases. 


FIG.  8:  VARIATION  OF  LOSS 
COEFFICIENT  WITH  BLADE 
NUMBER 
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ASYMPTOTIC  STABILITY  OF  PARTIALLY  DAMPED  MECHANICAL  SYSTEMS 
Gao  Weibin,  Beijing  Institute  of  Aeronautics  and  Astronautics 
Abstract : 

In  this  piper  the  asymptotic  stability  of  the  partially  damped  linear  me- 
chanicul  system  with  gyroscopic  forces  is  studied.  The  equation  of  motion  for 
the  system  is  a  second  order  linear  matrix  differential  equation.  The  stiffness 
matrix  is  assumed  to  be  positive  definite  representing  the  conservative  forces 
and  then  Lyapunov's  second  method  is  employed.  By  splitting  the  system  ma¬ 
trices  in  block  matrices  the  condition  of  asymptotic  stability  (e.  g,  the  con¬ 
dition  of  pervasiveness)  can  be  expressed  as  a  rank  condition  of  certain  ma¬ 
trix  with  dimension  lower  than  that  of  the  system.  Several  cases  with  different 
gyroscopic  coupling  are  discussed.  The  stability  condition  is  not  only  simpler 
than  that  obtained  by  former  authors,  but  also  gives  some  insight  to  the  me¬ 
chanical  properties  of  the  system  since  it  may  be  considered  as  two  coupled 
subsystems  then. 


I.  INTRODUCTION 

The  problem  of  stability  of  mechanical  systems  is  an  age-old 
problem.  In  recent  years  some  mechanical  system  problems  have 
emerged  in  aeronautic  and  astronautic  flight  instruments  and  other 
technologies,  again  arousing  the  attention  of  researchers.  For  the 
work  in  this  area  we  can  refer  to  literature  (1^-7)  . 

Taking  the  generalized  coordinate  vectors  of  mechanical  systems 

a3  C*i»  •“*  *•)»  then  the  equation  for  the  linearized  motion  dif¬ 
ferential  of  mechanical  systems: 

J43i+Dx+Gx+Kx-0  (1_1) 

is  a  second  order  linear  matrix  differential  equation.  In  it  M  ex- 

T 

presses  inertia  matrix,  generally  a  positive  definite:  M=M>0; 
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D  expresses  the  damping  matrix,  generally  a  positive  definite: 

D=DT>0;  G  expresses  the  gyroscopic  force  matrix,  it  is  a  slanted 
T 

symmetrical  G=G  ;  K  expresses  the  rigidity  matrix,  it  is  a  symmetrical 
T 

matrix:  K=K  .  Were  it  not  for  every  generalized  damping  force  not 
equal  to  zero  just  detD^O,  then  we  say  the  system  receives  incomplete 
damping,  or  the  system  is  called  partially  damped. 

/  \ 

In  literature  i  3)  with  the  classic  theorems  of  Kelvin,  Tait, 
Chetayev  ^  ' 

and  .  together  are  called  the  K-T-C  theorem:  the  stability  of 

system  (1-1)  is  equal  to  the  stability  of  the  conservative  mechanical 
system: 


Mx+Kx-  0 


(1-2) 


When  K)0(positive  definite)  and  E^O(semi-positive  definite), 
Zajac  points  out  if  at  this  moment  system  (1-1)  asymptotically  stab¬ 
ilizes,  then  the  damper  ie  called  "pervasive  damped".  Literature 
^4  and  5)  carry  out  research  on  this  condition.  To  the  gyroscopic 
system : 

tf*+0*+i£s- 0  (1_3) 

in  which  M>0,  n>0,  K>0,  literature  (S')  gives  the  criteria  for  asymp¬ 
totic  stability  (thus  conditions  of  pervasive  damping): 


order 


D(M-'K) 

DW'KY 


(1-4) 


and  literature  (7)  proves  the  criteria  for  asymptotic  stability  to 
even  more  generalized  mechanical  system  (1-1)  after  changing  into  a 
first  order  system: 


order 


CA  T  O  /.  -I 

h>.r* 


(1-5) 


The  results  of  this  text  (2-10)  and  (3-1)  are  somewhat  more  simp¬ 
lified  than  (1-5)  and  (1-4). 


% 

i 


II.  CRITERIA  FOR  ASYMPTOTIC  STABILITY 


K>0, 


To  research  mechanical  system  (1-1), 
T 

and  G=-G  ,  the  damping  matrix  can  be 


in  which  nxn  matrix  M>0, 
expressed  as: 


D 


detD.^O,  D, e#'". 


r<» 


(2-1) 


With  formula  (1-1)  written  into  a  first  order  form,  we  have: 

i-y,  i^-Af  tDy-Ar,Gy~M’lKx  (2-2) 

drawn  into  a  split  matrix  that  corresponds  with  (2-1): 


*  h  |,  > 


-(3  ^ 

ro,  oi  pA/,  A/,i  pG„  g„i 

A/-D-  ,  A/*'-  L  G- 

lD,  oJ  LA/,  A/J  Lg„  g4J 

rG,  G,1  {Kt  AT,1  fAT,,  /f„l 

A/-*G-J  L  A/*‘K-  L  *- 

LG,  Gj  U,  atJ  U„  a:J 

Thus  we  can  resolve  with  formula  (2-2),  to  obtain  two  systems: 


(2-3) 


G,yi-Gtyi-Ki*i-K*xt 
*,*y»  — Diy> — G,_y, — G ty i  -  K txt 

Constructing  Lyapunov's  function  to  system  (2-2): 


(2-4) 


It  clearly  is  a  positive  definite,  seeking  the  derivative  of  the 
solution  of  the  v  edge  (2-2),  we  obtain: 

«« — yTZ)y--yU>»yi  (2-5) 
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It  is  only  a  semi-degative  definite.  In  order  to  prove  that  (2-2) 
or  (2-4)  asymptotically  stabilise  and  seek  the  criteria  for  asymptotic 
stabilization,  we  use  La  Sail's  theorem  ^  .  For  this  reason  it  re¬ 
quires  satisfaction  of: 


teO  if  DQylS0 


(2-6) 


the  solution  of  formula  (2-4)  only  has  balanced  conditions,  x.j=0, 
x2=0,  y^O,  y2=0*  Condition  (2-6)  gives  out:  y^O  this  way  formula 
(2-4)  becomes: 


Oi  A,  Gtyt+Ktxt—  —KtA 
*! "^i,  X,"  —  Gty^  —  K sj4—  K txt 

In  which  A  is  a  constant  vector-product  constant, 
can  also  be  expressed  as: 


(2-7) 


Formula  (2-7) 


CAT,  j  *  —KtA 

LyJ 

H-r 0  7 

Ly.J  L  —  --C,Juv.J  L  A',J  Ly,J 

«T  0  ‘  1  ,-\°}a 

l -K<  .  -Gt  J  UJ 

With  formula  (2-8a)  we  seek  a  derivative  towards  t  to  obtain: 

C/C,  <?,3p **]-  0 

i»yj 

Substituting  formula  (2- 8b)  we  have: 


(2- 8a) 


(2- 8b) 


(2- 8c) 


(2- 9a) 


r* 


C/Ct  g.)  i 


Again  seeking  a  derivative  towards  t  to  obtain: 


(2- 9b) 


Continuing  down,  we  obtain: 

C/C, 


=•  -II:]' 


0,  Jk  “  0 ,  1,  2 » 


(2-9c) 


Matrix  L  is  2(n-r )x2(n-r )  dimensional.  According  to  the  Hamilton- 
Cayley  theorem  the  least  ^(n-r;  will  be  the  linear  combination 
I, *  * •L^(n~r)-1 ^  One  of  the  reasons  is  in  formula  (2-9c)  the  most  has 
2(n-r)  that  are  independent.  When: 

order  S“,2(«  — •’),  S* 


r  c/c,  (?tj 
C  Kt  GJL 


(2-10) 


l  C/c, 

formula  (2-9)  only  has  zero  solution:  *2=0,  y  =  0» 

y2~0»  x2=B 


(2-11) 


In  which  B  also  is  a  constant  vector.  At  this  time  equation  (2-7) 
becomes : 


k1a+k2b=o,  k3a+k4b=o 


Because  the  coefficient  matrix  of  the  above  equation  if  K  satisfies 
the  order,  one  of  the  reasons  is  A=0,  B=0.  This  precisely  shows  the 
system  only  has  zero  solutions  under  the  condition  of  (2-10),  thereby 
proving  formula  (2-10)  is  the  pervasive  condition  of  the  asymptotic 
stability  of  the  system. 


Now  we  prove  that  condition  (2-10)  is  also  necessary.  We  estab¬ 
lished  the  asymptotic  stability  of  the  system,  but  formula  (2-10) 
is  not  established.  At  this  time,  the  necessary  existing  vector  C* 
satisfies : 

C'EKtrS,  C*± pO  (2-12) 
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If  C*  satisfies: 


C/f,  GJL*C*  —  0 ,  A  -  0 ,  1,  2,  - 


(2-13) 


the  solution  of  the  differential  equation  (2-7)  possesses  the  form: 


*exp(Zi)Cf  L'H 


(2-14) 


Because  L  satisfies  the  order,  L~  exists.  According  to  literature 
(ll)  there  is 

-  J 

exp (£J)~  ^  a*(  /  )4* 

*•0 

In  which  a^.(t)  is  the  function  of  t.  With  a^(t)  multiplied  (2-13) 
all  formulas  later  mutually  increase,  to  obtain: 


H»C*-C/CuG,)exp(LOC*-  Q 


(2-15) 


Taking  integral  constant  C=C*,  substitute  formula  (2-15)  with  formula 
(2-14)  thus  to  obtain: 


“•  o  -C/C,  Gt)L  'l 


Again  substitute  form  (2- 8a)  with  (2- 8c)  and  consider  that: 


--n 


We  can  derive: 


K1A=G2K‘ 'k3a 


Because  satisfies  the  order,  and  the  right  end  is  not  equal  to 
K^A,  consequently  we  inevitably  have: 


A"  0.  /-«, 


(2-16) 
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This  clearly  shows  under  condition  (2-6),  formula  (2-7)  or  (2-8) 
has  no  zero  solution  (2-16),  and  is  contradictory  with  the  assumption, 
it  is  necessary  to  obtain  proof. 

The  asysmptotic  stability  fully  required  condition  of  the 
partially  damped  system.  Assume  the  system: 

M‘£+Djc+G£+Kx=0 

Where  M>0,  -G=G^,  K>0,  D=  |^°  ,  Do^0,the  order  of  the  rxr  matrix  Dt 

is  equal  to  r,  its  balanced  state  asymptotic  stability  fully  required 
condition  is  (2-10): 

order  S=2(n-r) 

Note  1.  Matrix  S  generally  has  2(n-r)  arrangement  and  in  liter¬ 
ature  7  matrix  (1-6)  has  2n  arrangement.  The  determination  of  the 
order  obtains  simplification,  this  particularly  is  notable  when  r 
approaches  near  n. 

Note  2.  Matrix  S  is  the  split  matrix  expression  of  each  coefficient 
matrix  through  system  (1-1),  one  of  the  reasons  is  that  it  is  easy  to 
analyze  the  influence  of  each  matrix  block  on  the  stability. 

Note  3*  On  direct  observation,  adding  the  fully  damped  Do  of  the 
system  that  is  expressed  by  x^  ,  inevitably  brings  about  x^(t)-*0,  and 
only  when  condition  (2-10)  is  established,  if  when  a  certain  gyroscopic 
coupling  and  rigid  coupling  exist,  then  can  it  cause  X£  to  also  have 
asymptotic  stability.  Condition  (2-10)  also  is  precisely  a  pervasive 
damped  condition. 

Note  4.  The  selection  through  the  coordinate  system  causes  M  to  be¬ 
come  the  diagonal  matrix,  and  will  cause  simplification  of  determination. 

III.  PARTIALLY  DAMPED  GYROSCOPIC  SYSTEMS  AND  N0N-GYR0SC0PIC  SYSTEMS 

1.  Non-gyroscopic  systems,  the  condition  of  G=0.  At  this  moment 
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the  system  equation  is: 


Mx+Dx+Kx-  0  (1-3) 

Considering  that  G=0,  condition  (2-10)  can  be  simplified  by: 


'  0 

/  ' 

CAT,  0)  ; 

rs.  01 

Kt  \ 

.-A', 

0  J 

S=  CO  K.J 

o 

y> 

\ _ 

.  S.A 

K.Kt 

j  —  CA‘.  A'4  O] 

K,K\ 

-CO  A'.  Aj' 

...  j 

I 

L  j 

In  which  certain  positive-negative  numbers  are  omitted  (this  does  not 
influence  the  determination  of  order),  to  obtain  equivalent  condi  'i 

order  So=n-r  (3 

Note  1.  Condition  (3-1)  has  comparatively  greater  simplifies 
than  condition  (1-4)  in  literature  M6).  This  is  similar  to  the  s. 
uation  in  the  second  section. 


Note  2.  We  still  can  directly  derive  condition  (3-1)  from  con¬ 
dition  (1-4).  By  formula  (2-1)  and  (2-3),  remember: 


D 


r  °i- 

lo  OJ 


O(AT’A')- 


'L\ 

.o 


Ll" 

*  ,  lA*D,Klt  L',*D,Kt 

OJ 


D(M'lK )*»  1  *  ,  Ll&L\Kt+L\K„  L\ f±L\Kt + L\Kt 

Lo  oJ 


sir 

L 


O  0 


} 


4?" ‘ ^ L%*K ,  +  L\‘lK r>  L?'±L\‘tKt+LrtK. 


They  have  the  following  equal  relationship: 
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Order  S 

D,  0 

0  0 

= 

-  D,  0  I 
L\  L\ 

=  order 

On 
0  L\ 

=  order 

L\~L\ 

0  0 

order 

.  o  ir- 

ir  Lyi 
o  o 


=r+order 


r  «  i 

r  d,k,  •> 

rDtK,  h 

Li 

”  r  +  order  j 

I  lik4 

—  r  +  order 

D,KtK  4 

-  L;-xKt  . 

r  +  order  So 


The  3.5.  and  6  equal  formulas  of  the  above  formulas  go  through  multiple 
determinant  transformations  to  be  obtained,  the  most  unequal  formula 
transformation  is  lengthy  and  the  process  is  omitted.  This  way  it 
proves:  "Order  S=n"  is  equivalent  to  "Order  So=n- r" . 


Note  3.  According  to  the  dissipation  we  can  conduct  the  following 
classifications  with  mechanical  systems: 

a.  Fully  damped,  D>0;  System  asymptotic  stabilization. 

b.  Partially  damped,  D^O: 

(1)  Pervasively  damped,  order  So=n-r:  System  asymptotic  stab¬ 
ilization  . 

(2)  Not  pervasively  damped,  order  So<n-r:  System  stable  but 
not  asymptotically  stable. 


2. 


Partially  gyroscopic  system: 


O' 

f 

o 


the  dimension 


of  block  matrix  Go  is  equivalent  to  the  dimension  of  Do.  In  this 
condition,  M_1G=|G1  J  ,  condition  (2-10)  also  becomes  (3-1),  G, 

\  3  °3 

does  not  produce  an  effect  on  the  stability.  This  can  be  expected, 
because  the  system  that  x^  expresses  is  fully  damped. 


3. 


Partially  gyroscopic  system: 


In  this 


condition,  due  to  inertia  coupling,  adding  gyroscopic  force  of  the 
system  that  expresses  will  produce  effect  on  the  whole  system.  It 


1  04 


m 


influences  the  order  of  the  S  matrix,  if  influencing  pervasive  damped 

T 

conditions.  But  when  inertia  coupling  doe3  not  exist,  M^M^O,  then 
condition  (2-10)  becomes  condition  (3-1 K  and  proof  is  omitted. 


4.  Gyroscopic  coupling  system: 


This  moment  generally  speaking,  the  coupling  of  gyroscopic  forces 
greatly  influences  the  pervasiveness  of  damping,  and  moreover,  even  if 
the  inertia  force  is  without  coupling,  by: 


'G, 


G.'l  r  O  Afl'G'.'l 

OJI-MVGI  0  J 


We  can  obtain: 

fc 

It  is  thus  evident  that  <?,  coupling  has  influence  on  the  pervasiveness 
of  the  damping  of  the  system.  It  can  cause  partial  damping  of  the 
system  to  become  pervasive.  From  this -we  can  make  the  system  asymptot¬ 
ically  stable. 


IV.  EXAMPLES 


1.  Consider  the  mechanical  system  that  is  shown  in  Fig.  1. 


*  ffft  **■  S  4  **■!  ,  %•!  . 


FIG.  1:  MECHANICAL  SYSTEM  WITH 
PERVASIVE  DAMPING 


The  equation  for  the  system  is: 


M'x+Dx+Kx -  0 
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diagCm,,  ••*, 

A.  +  A.  -A,  _  0 

1  d,  +  dt  —dt  _  0  ! 

—  k,  At+A, 

-d,  d.  +  dj 

0 

• 

1 

-«? 

+ 

■ 

.  D=  —  d».t 

•  *. 

1 

!  0  -d.-,  dm-t  . 

-=? 

-a 

1 

O 

i  0 

CO 

CO 

According  to  the  definition,  we  have: 


fA\ 

A:' 

0 

0  1 

(  A„. ,  4-  A. 

_ k,  '■ 

| 

II 

£ 

* 

* 

"ll 

l 

m.., 

m.., 

u. 

a4. 

i  0 

0 

9 

H  *. 

t 

J. 

I 

1  — 

0  i 

1  m. 

m.  ; 

J 


Due  to  : 


r  K*  ^ 
■ft1  A'X  : 

L  ...  J 


o 


2 


We  can  know  the  system  has  asymptotic  stability,  damping  is  pervasive, 
and  the  conclusive  physics  phenomenon  is  obvious. 

2.  Consider  the  mechanical  system  that  Fig.  2  shows. 


FIG.  2:  MECHANICAL  SYSTEM  WITH 
NONPERVASIVE  DAMPING 


The  equation  for  the  system  is: 


MX+D*+Kx-  0 


At  +  *| 

Af“diag(m,  •••  m4),  K  ”  — 

D“diag(0  dt  0  0),  q 


-k,  O  O' 
A.+A,  -A,  0 

—  A,  A,+A,  A, 


O  O  -*«  A* ) 
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Changing  the  sequence  of  coordinates  and  x the  coefficient  matrix 
separately  becomes: 


M  —  diag(m,  m,  m,  m4),  K  — 
D-di«g(d,  OOO), 


k,+k,  —  fc,  ~l*t  O 

—  ht  hl  +  kt  O  O 

—  k,  0  +  —  /*4 

O  0  -A«  A. 


Computations  are: 


.K, 


h 1 4*  kj 

m. 

0 


O 


o 


k,  -I-  kt 

m , 

_ k±_ 

~  m. 


O 

k, 

m, 

A_ 

1914 


Order  So=order 


(  _A_ 

m. 

k3 

m. 

> 

0 

1 

A_ 

k,  _ 

4*  kA 

Jh  kt_ 

m. 

m. 

m. 

mz 

m.  ms 

_/A±V 

f  k'~ 

_ji.( 

k,+kt  y 

fc-fe4  fe4 

\  rn,  y 

mt 

m.  \ 
kak 

m3  ) 

4 

m.m,  m, 

.  '  k,ki 

If  k^=k2=k^=k,  m^  =m2=ni^=m^=m ,  then  S  descending  order  condition  is 
k=3/2  k^.  At  this  moment*  partial  damping  is  particularly  pervasive, 
and  possible  motion  emerges:  m2  is  unmoving,  the  motion  of  m^  with  the 
motion  of  m^  and  m^  are  reversed  resonant  motion. 


V.  CONCLUDING  REMARKS 


With  classic  mechanical  system  differences,  there  the  primary 
research  is  on  conservative  systems,  and  in  modern  times  in  mechanical 
systems  artificial  dampers  are  often  recommended  and  are  regarded 
as  the  composed  means.  We  recommend  dampers  of  comparatively  smaller 
numbers  in  order  to  ensure  that  the  problem  of  system  asymptotic 
stability  has  certain  reference  significance.  When  analyzing  partially 
damped  mechanical  systems,  with  the  coefficient  matrix  of  the  system 
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a  split  block,  we  can  establish  an  ordered  criteria  of  a  lower  di¬ 
mension,  and  also  it  is  convenient  to  analyze  the  mechanical  properties 
of  the  system. 
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AERODYNAMIC  COEFFICIENT  IDENTIFICATION  OF  TIME-VARYING  AIRCRAFT 
SYSTEM  AND  ITS  APPLICATION 


Wang  Tong,  Beijing  Research  Institute,  China  Precision  Machinery 
Corporation 

Abstract : 

Aerodynamic  coefficient  identification  for  time-varying  aircraft  system  has 
been  studied.  Usually,  the  time-varying  system  identification  is  quite  difficult. 

On  the  basis  of  practical  measurement  in  flight  tests  and  by  analyzing  trend  of 
the  coefficients  in  aircraft  time-varying  differential  equations,  it  is  possible  to 
transform  the  individual  time-varying  coefficient  into  a  known  function  multi¬ 
plied  by  an  unknown  constant.  These  unknown  constants  are  referred  to  as  unde¬ 
fined  coefficients.  With  the  aid  of  the  Newton-Raphson  method  extended  to  the 
time-varying  coefficient  differential  equations  the  undefined  coefficients  can  be 
evaluated  by  iteration  calculation.  In  this  way  the  complicated  time— varying 
aircraft  identification  can  be  carried  out. 

In  order  to  verify  this  calculation  the  aerodynamic  coefficients  of  an  aircraft 
have  been  evaluated  by  the  data  taken  from  one  of  its  unsteady  flight,  and 
the  reliability  of  aerodynamic  coefficients  obtained  from  this  identification  has 
been  discussed.  The  calculation  results  matched  very  well  the  test  data  and  made 
a  contribution  to  the  improvement  of  aircraft  fight  tests. 

I.  PREFACE 

According  to  the  identification  theory,  using  flight  test  data 
to  compute  the  aerodynamic  coefficient  of  the  flight  instrument  has 
a  very  practical  significance.  In  order  to  use  the  identification 
theory  to  obtain  a  more  accurate  aerodynamic  coefficient,  except  for 
the  problwm  of  computation,  we  still  must  provide  a  series  of  ad¬ 
vantageous  conditions:  for  example  the  aircraft  has  suitable  motorized 
motion,  causing  the  flight  test  data  to  include  sufficient  aerodynamic 
coefficient  information;  it  requires  the  arrangement  of  suitable  more 
precise  measurement  bearing;  the  aircraft  mathematical  model  is  fairly 
close  to  realistic  conditions.  For  this  reason  we  can  carry  out 
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specialized  flight  tests,  after  the  aircraft  flight  stabilizes,  the 
rudder  surface  carries  out  movement  according  to  a  pre-determined 
course,  the  aircraft  makes  motorized  flight,  bearing  the  corresponding 
data  that  has  been  recorded.  Under  this  kind  of  condition,  from  the 
flight  test  data  we  compute  the  aerodynamic  coefficient  of  the  air¬ 
craft,  the  motion  process  of  the  aircraft  generally  is  a  constant 
coefficient  differential  equation. 

Real  flight  tests  sometimes  meet  with  problems  and  unstable 
flight  occurs.  If  the  problem  does  not  stem  form  the  pilot  apparatus 
or  other  components,  then  it  is  possible  that  the  original  mathematical 
model  of  the  aircraft  is  not  accurate.  The  originally  determined 
aerodynamic  coefficient  and  the  real  aerodynamic  coefficient  in  flight 
are  not  in  accord.  This  was  objectively  then  it  is  necessary  to 
utilize  flight  test  data  t-o  compute  the  aerodynamic  coefficient  of  the 
aircraft.  Clearly  at  this  time  the  conditions  for  computation  are 
insufficient,  the  rudder  motion  was  not  specially  designed  for  the 
purpose  of  obtaining  the  aerodynamic  coefficient,  and  there  is  problem 
motion.  The  measured  parameter  that  is  advantageous  to  parameter  id¬ 
entification  has  no  certain  provision  to  compute  the  aerodynamic  co¬ 
efficient,  we  can  only  use  the  existing  measured  parameters  of  flight 
tests,  particularly  the  linking  together  of  frequent  problems  and 
acute  transformation  of  parameters.  Therefore  we  cannot  take  the 
motion  equation  of  the  aircraft  and  treat  it  as  a  constant  coefficient 
differential  equation,  but  should  treat  it  as  a  variable  coefficient 
differential  equation.  This  text  researches  parameter  identification 
under  this  kind  of  condition. 

II.  EQUATION  FOR  VARIABLE  COEFFICIENT  AIRCRAFT  SYSTEM 

The  following  only  researches  the  side  motion  equation  of  the  air¬ 
craft,  the  side  motion  compared  to  the  vertical  motion  is  a  little 
more  complicated,  and  hte  possibility  of  problems  that  emerge  is  some¬ 
times  a  little  more.  Additionally,  consider  that  the  aircraft  is  an 
elastic  body,  due  to  elastic  vibration  of  the  aircraft,  in  flight  tests 
the  angle  degree  and  angle  velocity  signal  that  the  telemetering 
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transducer  is  affected  by,  except  for  including  the  angle  degree 
and  angle  velocity  signal  of  aircraft  real  flights,  still  includes 
the  rising  angle  degree  and  angle  velocity  information  of  the  air¬ 
craft  elastic  vibration  in  the  transducer. 

The  equations  for  aircraft  motion  are  as  follows: 

*  )*i  +  -4u(  t  )*,  +  .!„(  t  +  t  >,  +  £„(  t  )u,  'i 
=  *  )*i  +  -4;!(  t  t  )x3 t5.,(  t  )ul  +  Btt(  <  )“» 

a  (  t  ) 

+  *t  +  A33(  t  )X,  +  AU(  t  )x, 

*«-*!  (1) 

*»**» 

+  t  )xt  +  D3(  t  )*,  +  D,(  t  )ut  ' 

*,*■*,  +  *, 

*,-*,+*,  ) 

Where 

is  trundle  angle  velocity 
X£  is  yaw  angle  velocity 
x^  is  side  slip  angle 
x^  is  trundle  angle 
x^  is  yaw  angle 

X/  is  signal  of  the  elastic  vibration  of  the  aircraft  on  the  yaw  angle 
"  measurement  (rudder  free  gyration) 

Xj  is  signal  of  the  elastic  vibration  of  the  aircraft  on  yaw  angle 
'  velocity  measurement  (rudder  damped  gyration) 

Xg  is  signal  of  included  aircraft  elastic  vibration  on  rudder  damped 
0  gyration 

Xg  is  signal  of  included  aircraft  elastic  vibration  on  rudder  free 
gyration 

u.|  is  aileron  deflector  angle 
Uj  is  yaw  rudder  deflection  angle 

A.j  i  (t) ,  A^j(^)*  A^j' t) ,  A  ( t ) ,  A22^)>  A  23  ( t ) ,  A  ^  ^  ( t ) ,  A3  ^ ( t )  are 
the  aerodynamic-  coefficients  of  the  aircraft 

B^(t),  B2i(t),  are  rudder  surface  aerodynamic  coeffic¬ 

ients  of  the  aircraft 


^  is  aircraft  elastic  vibration  damping  coefficient 
a>  is  aircraft  elastic  vibration  frequency 
a(t)  is  angle  of  attack 

D2(t),  D^(t),  D^(t)  separately  are  yaw  angle  velocity,  side  slip  angle, 
yaw  rudder  angle  regarding  the  elastic  vibration  coecion  force  of  the 
aircraft. 

Computing  aircraft  aerodynamic  coefficient  from  flight  test  data 
(if  we  seek  Aij(t),  Bij(t),  and  Di(t)),  because  these  parameters  all 
follow  time  variations,  moreover  the  amount  is  fairly  much,  therefore 
they  are  difficult  to  directly  obtain.  The  following  makes  suitable 
transformations  to  these  parameters,  with  appropriate  identification 
computations . 

Coefficients  Aij(t),  Bij(t),  and  Di(t)  receive  the  influence  of 
many  factors.  If  altitude  H,  flight  Mach  number  M,  and  aircraft  weight 
G  or  inertia  I  that  rise  due  to  fuel  consumption  transform,  then  we 
have  the  expressed  formula: 

A,,i  i  )-v4,,U/(  O.  ),  Gi  t  )) 

B„i  t  )-£,,CAf(  t  ),  t  ),  Gi  t  )) 

D,(  t)-D,CMi  t  ),  Hi  O,  G CO) 

But  these  factors  M(t),  B(t),  and  G(t)  in  flight  tests  are  com¬ 
paratively  easy  to  measure  or  to  obtain  by  computation,  and  are  com¬ 
parable  with  real  flight  aerodynamic  coefficients  that  are  difficult  to 
ascertain.  The  believeable  estimate  of  the  value  of  these  factors  must 
be  very  large,  therefore  in  the  process  of  flight  test  data  that  is  used 
to  compute  the  aircraft  aerodynamic  coefficient,  we  can  extract  these 
already  known  time  varying  factors,  and  the  remaining  unknown  sections, 
taking  undefined  coefficient  expression,  in  this  way  Aij(t),  Bij(t),  and 
Di(t)  can  be  written  as  Aij (t)=aijFij (t) ,  Bij (t)=bijGij (t) ,  Di(t)=diHi(t) . 
In  which  ai j ,  bij,  and  di  are  undefined  coefficients,  F^jCt),  Gij(t), 
and  Hi(t)  are  already  known  time  varying  functions  with  computed  A^Ct)* 
=a^F11(t)  as  an  example  to  explain  this  point. 
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Where 

q  is  velocity  head  q=q^M(t),  H(t)^; 
v  is  velocity  V=VfM(t),  H(t)^  ; 

is  aircraft's  largest  cross  section; 
is  aircraft's  largest  diameter; 

Jx  is  elastic  axis  directional  rotation  inertia  Jx=Jx(t); 
m^x  is  damping  coefficient  of  the  rolling  force  moment,  by  theoretical 
computaiton  or  wind  tunnel,  it  follows  the  regulation  of  the  flight 
Mach  number  M  is  already  known,  therefore  we  have  m^x=a^f(M(t)^  ,  a^ 
is  an  undefined  coefficient,  f(M(t^|  is  the  similar  formula  of  m^x  fol¬ 
lowing  M  variable  regulation. 


Therefore : 


A,(  * ) 


5 1  CAf(  /  ).  H(t  )35tDl 
”  t  )FCA/(  i  ),  //(  t  )1 


a,,/CA/(  I  '  ) 


Using  the  same  kind  of  method  to  obtain  the  rest  Aij ( t) =aijFij ( t) , 
Bij (t)-bijGij (t),  Di(t)=diHi(t) . 


Taking  the  above  formulas  and  substituting  equation  (1) 


we  thus  have: 


t  )*1+a„F11(  <  )jt|+a,jFu(  t  )x,+b,,G,,(  <  )«,+6|,(?t,(  i  )#,"! 

I  )*i  "bOu^uf  t  t  )*,  +  <*„<?„(  I  )u, +6..Cr,s(  t  )#, 

a  (  (  ) 

*«  ”  57  3  -*!+■*,+ Oj,F „(  t  )x,+  auFu(  I  )xt  < 

*4-*.  1 

r,  =  X7  | 

x,=  +  «/,//.(  t  )*,  +  <*,//,(  /  )x%  +  dxH£  t  )tt, 

-r,  =  x,  +  x, 

*,  ~  x,  t*. 


(2) 


Therefore  we  compute  the  aerodynamic  coefficient  from  flight  test 
data  and  the  problem  of  i  and  a  in  a  nutshell  is  tc  seek  the  undefined 
coefficients  aij,  bij,  di,  i  and  ® of  the  variable  coefficient  differen- 
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tial  equation  (2). 

III.  COMPUTATION  FOR  PARAMETER  IDENTIFICATION 

The  above  already  establishes  the  mathematical  model  for  the 
variable  coefficient  aircraft  system.  In  generalized  aircraft  tests, 
through  a  telemeter  we  can  obtain  ^i,  u,,i,  x3ineasured,  x9imeasured. 
and  x,i  ,  (the  lower  "i"  symbol  is  a  time-order  of  telemetered 

sampling).  The  task  of  parameter  identification  is  simply  to  take  the 
telemetered  rudder  yaw  u^ i  and  U2i  as  input,  select  suitable  aerodynam¬ 
ic  coefficients,  substitute  these  aerodynamic  coefficients  in  equation 
(2),  solve  the  square  and  as  the  smallest  of  the  difference  of  x^i, 

Xgi,  Xgi,  and  x^i  tha  equation  (2)  caused  to  obtain  and  the  corres¬ 
ponding  telemetered  values  x.i  „  ,,  x0i  , ,  xni  and 

^  6  1  measured  8  measured  9  measured 

x41measured‘ 

In  the  mathematical  expression,  if  we  establish  a  function  J 


(measured). 


(3) 


take  the  undefined  coefficients  aij,  bi j ,  di,  t,  and  a>  ,  and  compose 
a  vector  c. 

cT*  C  *u,  ®u»  «ii,  6,,,  6,„  o„,  a,„  a,„  6tt,  a,„  o,4,  fc,  ®»  d „ 

dit  d,J 

(T  expresses  the  turning  position). 

J  value  is  the  function  of  vector  c,  by  J  value  towards  vector  c 
we  seek  the  smallest  value  to  determine  vector  c. 


By  using  the  Newton-Raphson  method  of  iteration  formula  we  can 
obtain  the  J  v  ,?.ue  to  reach  the  smallest  vector  c  ^  . 


114 


Ciui 


(4) 


Where  Xi  is  the  vector; 

Xi“<x1i»  x8i ’  V*  x4i} 

T 

X  =ix  i  xi  xi  xi 

imeasured  v  1  measured'  8  measured'  9  measured'  4  measured) 

X7cXi  expresses  the  differential  of  vector  Xi  to  vector  c. 


By  formula  (4)  we  can  know,  the  iteration  computation  of  vector  c 
can  in  a  nutshell  be  computed  VcXi.  Because  the  lower  symbols  that  the 
computations  used  to  be  reached  are  comparatively  many,  for  clear  ex¬ 
pression,  the  following  computations  omit  the  lower  synbol  "i" . 


VcX  uses  the  matrix  expressed  as: 


V.*- 


i 

dan 

*(Xt  +  X7) 
da,  | 

*(*»  +  *.) 

**« 

.  *an 


*>i, 

*(xt  +  Xy) 

*an 

<>(*.  +  *.) 

datl 

_jxi_ 

**u 


■ 

3a„ 

*(*,  +  *?) 

dais 

*(*.  +  *.) 

*®i» 

*xt 


to,  to,  **i 

dbn  *btt  **t, 

<(x,  +  xT)  *(*»  +  *t) 

^6,,  ^®ti 

<>(*«+*«)  *(x.+xt)  a(x,+*,) 

06,,  **11 

to,  to,  to, 

a6„  *6,»  *°«i 


to, 

to, 

St" 

to. 

*»» 

to,, 

*6 

St 

ia,, 

to,. 

*(*«  +  X,) 

«l« 

*(X,  +  X7) 
to,, 

*(x,  +  X7) 
dbtl 

d(x,  +  x7) 

*(X,  +  XT) 
to,, 

*(*»+*») 

tosi 

*(*.+*< 

«n 

iL  . 

<»(xt+xt) 

to,, 

*(*»+*»>  . 
*6„ 

<»(x,  +  x,) 
d6„ 

<»(x,+x,) 

to,. 

<Kx,  +  x,) 
to,. 

to. 

to. 

to. 

to. 

«M 

to,. 

*b» 

to,. 

to,. 

0 

0 

0 

0 

0 

> 

-?r 

toT 

to) 

dx7 

»dx 

to, 

**» 

■  A. 

H 

to, 

to> 

to, 

»dt 

O' 

0 

0 

0 

9 

i 

(5) 
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Each  element  of  matrix  ^cX  is  obtained  by  resolving  the  differential 
equation. 


The  method  of 


ax, 


ax. 


»x. 


ax. 


ax. 


ax. 


'  *°H  *  *®!i  '  dOu  '  dOu  *  *»W 

Resolving  the  following  differential  equations: 

-%r- •••F" ( '  Hs :-fr+°“F" < 1  >Str+‘"F»  < 1  >-£*- 


+  /ru(  t  )x 


6*± - c  '  '  dX-l—j.»  c  (  #  **\ 


aa, 


t  y 


da, 


+  otlFit(  t  )- 


da,, 


+  o„Ft,  (  / 


da, 


dx3 


da„ 

57.3 

dX4  _ 

d*.__ 

da,, 

da„ 

_d*»  _ 

dx. 

da„ 

da„ 

**.  . 

dX, 

da,, 

do  u 

da, 


da, 


da, 


da, 


dx. 


-*-2Zr-*-2r+*#t  <  *  )-» 


dX, 


dan 


dan 


da. 


The  coefficients  of  the  above  equation  «n>  a„4  a„,  a„,  aM>  a„,  a,4> 

4,  .  d..d,  are  the  elements  of  vector  c.  In  this  iteration  method  we 

take  the  value  of  c^  to  go  through  computations  to  obtain  the  value  of 
°K+1 ’  therefore  in  resolving  the  above  equations,  these  coefficients 
are  already  known  values. 


When  assuming  the  prior  conditions  t-t  i 


„  JZfVjl-  »  dX,(j,)  =  **  ,({.)_  _  **,(/,)  _  dX,(t,)  _  dX.ff,)  _ 

da,,  da,.  da. .  da,,  ~  ~  “  ”  -  "  ® 


da. 


da, 


If  when  t=tQ  flight  trajectory  is  even  and  steady  and  variation  is 
slow,  this  kind  of  assumption  is  permitted. 

Resolving  the  above  equation,  we  then  can  obtain: 


dX, 


dx. 


dx. 


dX, 


dx. 


dX , 


d°i i  '  da, ,  da,,  '  do,,  '  da,,  '  da,, 
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Using  a  similar  method  we  can  obtain  the  yaw  derivative  of 

V  *2-  V  V  V  an(i  x7  t°"ar'13  a12*  a13‘  b11  ’  b12-  a21  ’  a22* 

a23’  ^21*  ^22’  a33*  anc*  a34’  anc*  ^aw  derivative  °f  and  x^  to¬ 
wards  |  ,  co  ,  d^  ,  and  d^. 

After  obtaining  VcX,  substitute  equation  (4),  we  then  can  use 
the  iteration  method  to  obtain  vector  c. 

IV.  APPLICATION 

The  following  practical  example  from  a  flight  test  explains  the 
application  of  the  above  variable  coefficient  aircraft  system  identif¬ 
ication  theory  and  formulas.  Certain  aircraft  in  the  operation  period 
emerge  form  steady  flight  and  develop  unsteady  flight,  so  much  that  the 
flight  diverges  (See  Fig.  1  through  Fig.  5).  Through  investigating, 
the  factors  of  unsteady  flight  are  located  in  the  aircraft's  original 
mathematical  model  that  is  not  accurate,  if  the  aerodynamic  data  in  the 
wind  tunnel  and  the  real  aerodynamic  data  of  the  aircraft  in  flight 
test  do  not  correspond,  therefore  it  is  necessary  to  compute  the  air¬ 
craft  aerodynamic  coefficient  from  the  flight  test  data,  and  find  the 
existing  problem  of  the  flight's  instability. 

The  above  aircraft  mathematical  model  is  a  variable  coefficient 
differential  equation  that  includes  elastic  vibration  (See  equation  2). 

The  numeric  value  of  the  realistic  aircraft  system  substituting 
with  equation  (2),  the  aerodynamic  coefficients  that  are  necessary  to 
compute  are  regarded  as  undefined  values,  and  we  can  obtain: 
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i,-a,,(0.26  t  -O.lTK+oJO.’e  t  -0. 17)*,+ a,  ,(0.46  t  -rl.06)x, 

+  bu( 0.46  f  -1.06)*,  +  6,.(0.46  t  -1.06)u, 

*i“Oii(0.26  t  -0.17)*,  +  a.,(0.26  I  -0. 17)x,+a„(0.46  t  -1.06)  *, 

+  bu  (0.46  t  —  1.06)u,  +  6„(0.46  t  -1.06)tt, 

*»- <3^--«l  +  *»+a»t(0.26  t  -0.17)*,+  *, 

*«=*T 

x,--|<0x,-«>*xl  +  </1(o.4  /  -0.82)at  +  rft(0.4  <  -0.82)*, 

+  of, (0.4  /  -0.82)*, 

*,—*,+*, 

*,=»*,+*,  J 

From  flight  test  telemetered  data  we  can  obtain  rudder  yaw  and 
u^.  In  equation  (6),  u^  and  u^  are  regarded  as  input  symbols,  therefore 
in  identification  computation,  it  is  an  already  known  value. 

Flight  test  telemetered  data  also  gives  angle  degree  and  angle 
velocity  sybols  x, measured,  x8measured>  x9measursd,  x^easured. 

Applying  formula  (4).  seeking  the  undefined  coefficient  in  for¬ 
mula  (6),  causes  the  square  and  J  as  the  smallest  of  the  difference  of 
,  Xg,  Xg,  and  x^  that  are  computed  in  formula  (6)  and  the  measured 

values  x-j measured ’  x8measured’  x9measured '  an<^  x4measured* 

The  computations  in  the  computer  clearly  express  that  the  undef¬ 
ined  coefficients  are  many,  they  not  only  greatly  increase  computer 
time  but  also  sometimes  in  the  iteration  process  it  is  not  easy  to 
converge.  The  preliminary  computations  of  the  example  in  this  text 
indicate  that  the  above-mentioned  aircraft  elastic  vibration  only  holds 
a  very  small  amount  in  the  exhaust  numeric  value  of  the  yaw  angle  sen¬ 
sor  and  yaw  angle  velocity  sensor.  Under  the  premise  that  does  not 
influence  attainment  of  the  aerodynamic  coefficient  value  that  plays 
a  major  role  towards  the  instability  process,  the  computer  did  not  take 
4  »  ©  d.| ,  d^t  dy  and  regard  them  as  undefined  values,  the  damping 
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coefficient  £  and  vibration  frequency  •  take  the  test  value,  and  , 
d£t  and  d^  take  the  designed  value.  By  the  aerodynamic  coefficient 
that  the  computations  obtained  we  see  the  following: 
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Taking  the  aerodynamic  coefficient  that  the  computer  obtained  and 
substituting  them  with  equation  (6),  we  take  the  yaw  u^  and  u^  that  the 
telemeter  obtained  as  input,  resolve  equation  (6),  and  computed  value 
and  flight  test  value  can  better  agree  (See  Fig.  3  through  Fig.  5). 

When  practically  applying  these  aerodynamic  parameters  that  id¬ 
entification  computations  obtained,  because  telemetered  parameters 
always  have  errors,  the  bea-ring  of  this  example  of  flight  test  did  not 
in  advance  consider  the  requirements  of  identification  computations, 
and  errors  were  also  somewhat  large.  Clearly  we  only  must  have  those 
aerodynamic  coefficients  that  were  not  very  sensitive  to  the  telemeter¬ 
ed  signal  errors,  then  they  can  be  believable.  If  these  aerodynamic 
coefficients  are  closely  linked  together  with  the  appearance  of  the 
aircraft  of  the  whole  flight  test  trajectory  transformation,  it  has 
a  significant  influence  on  the  flight  test  results,  therefore  even  if 
the  telemetered  data  has  some  error,  after  identification  computations, 
these  aerodynamic  coefficients  also  cannot  have  a  great  many  transfor¬ 
mations  . 

In  order  to  determine  which  aerodynamic  coefficients  possess  re¬ 
lative  stability  toward  the  transformation  of  telemetered  data,  we 
made  the  following  investigation:  Taking  the  aileron  deflation  uimeasureci 
and  the  yaw  rudder  deflection  u2measured  example  as  the  input 

signal,  with  each  different  disturbance  on  x1ffleasured,  x8aeasured. 

xn  j,  and  x,  ,  later  take  the  telemetered  value  of  this 

^measured  4measured 

kind  of  artificial  disturbance  and  regard  it  as  a  supposition  measured 
value,  and  again  form  these  supposition  values  carry  out  aerodynamic 


coefficient  identification  computations.  The  computer  results  show 
the  rising  trajectory  that  corresponds  to  the  rudder  deflection  trans¬ 
formation  in  this  example,  after  each  different  disturbance,  some 
aerodynamic  coefficient  variations  are  fairly  large,  but  the  aileron 
aerodynamic  coefficient  b^  and  yaw  rudder  aerodynamic  coefficient  b22 
however  have  corresponding  stability.  This  explains  the  value  of  b^ 
and  b  .^possess  certain  believable  degrees. 


The  rudder  surface  aerodynamic  coefficients  that  identification 
computations  obtain  are  much  greater  than  originally  designed  values, 
causing  amplification  system  of  the  aircraft  and  pilot  apparatus  return 
loop  to  greatly  raise  the  original  designed  value.  This  then  explains 
the  factor  of  aircraft  instability.  In  later  test  flights,  we  used 
reduced  pilot  apparatus  return  loop  to  amplify  the  coefficient  in  an 
equal  manner,  causing  the  aircraft  to  obtain  a  successful  stable  flight. 


V.  CONCLUDING  REMARKS 


This  text  researched  the  problem  of  aerodynamic  coefficient  iden¬ 
tification  of  variable  coefficient  flight  systems  with  elastic  effect. 
Through  the  analysis  of  variable  coefficient  transformation  regulations 
on  inside  equations,  it  turned  the  variable  coefficient  in  the  equation 
into  a  variable  regulation  already  known  function  and  unknown  constant 
product,  regarding  these  unknown  constants  as  undefined  coefficients, 
applying  the  Newton-Raphson  method  to  reach  the  variable  coefficient 
equation,  and  through  iteration  computation  seek  these  undefined  co¬ 
efficients.  This  kind  of  method  this  text  used  is  generally  so  that 
the  parameter  identification  computation  of  the  corresponding  difficult 
variable  coefficient  aircraft  system  can  be  achieved.  This  text  applies 
the  above  method  to  successfully  compute  its  aerodynamic  coefficient 
from  the  flight  test  data  of  one  unsteady  flight  of  an  aircraft.  When 
applying  in  practice  the  aerodynamic  coefficients  that  the  identifi¬ 
cation  computations  obtained,  if  the  measurement  errors  are  large,  we 
must  pay  attention  to  the  believable  degree  of  the  aerodynamic  co¬ 
efficient  that  these  computations  obtained  in  analysis  and  differen- 
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tation. 


This  work  was  carried  out  under  the  guidance  of  Comrade  Bao 
Kerning.  In  the  process  we  obtained  much  help  from  Comrade  Wang 
Guolin.  Please  accept  my  sincere  thanks. 


FIG.  1:  THE  TIME  HISTORY  OF  FIG.  2:  THE  TIKE  HISTORY  OF 
AILERON  DEFLECTION  MEASURED  RUDDER  DEFLECTION  MEASURED 
IN  FLIGHT  IN  FLIGHT 


IN  FLIGHT  AND  COMPUTED  IN  FLIGHT  AND  COMPUTED 
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AN  ADJUSTMENT  METHOD  AND  ENGINEERING  REALIZATION  FOR  CONTROL 
CURVES  OF  THE  TWO- VARIA3LE  FUNCTION 


Shao  Rongshi 

Design  Department,  Shenyuan  Aircraft  Corporation 
Abstract : 


Thmugh  designing  a  bypass-door  system  for  a  supersonic  aircraft,  which 
implemented  adjustment  in  conformity  to  the  control  curves  of  the  two-variable 
function  (Fig.  1  ),  a  method  for  engineering  realization  of  some  separable  mul¬ 
tivariable  functions  has  been  put  forward,  It  means,  for  any  control  function 
of  n  variables  /  (jc„  x„  •••,  jc.),  if 

-Zli* ir  *».  =  n 


dXidX, 


j  =  1  ,  2  ,  • 


the  function,  called  separable  by  the  author,  can  be  divided  into  subfunction  of 
each  single  variable,  then  the  sum  of  these  subfunctions  is  just  equal  to  the  ori¬ 
ginal  function,  i.  e. 

/(* i,  *i,  — ,  x.)“/i(*i)  +  /t(*t)  +  —  +  /•(*•). . 

A  practical  bridge  circuit  with  parameters  realizing  the  control  curves  for  the 
bypass-door  has  been  designed  (see  Fig.  4). 

Moreover,  the  simple  method  of  changing  the  characteristics  of  the  control 
curves  by  adjusting  resistors  on  the  bridge  arms  has  been  adopted  according  to 
the  test  flight  requirements  of  the  svstem,  and  an  actual  numerical  example  has 
been  given. 

The  above  mentioned  engineering  realization  and  adjustment  method  for  the 
two-variable  function  have  done  their  bit  for  the  research  on  the  simulator 
reproducing  the  function  of  more  than  two  variables  and  may  be  also  valuable 
for  design  of  the  system  with  multiple  input  but  single  out-put. 


I.  INTRODUCTION 

The  design  of  many  automatic  adjustment  systems  in  modern  aircraft 
design  requires  a  controlled  target  according  to  predetermined  multi- 
variable  function  curves  to  carry  out  open  formula  adjustments.  For 
example,  the  air  intake  passageway  of  the  high  speed  fighter  plane 
requires  according  to  M  number  and  flight  angle  of  attack  these  two 
parameters  to  adjust  the  air  intake  cone or  oblique  plate;  the 
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motorized  wing  flap  of  the  whole  front  edge  of  the  wing  spread  of 
America's  light  duty  combat  plane  YF-16  follows  angle  of  attack  in¬ 
creases  and  automatic  lower  incline,  and  follows  the  increases  and 

(2^ 

automatic  upper  convergence  of  the  M  number^  1  .  The  key  problem  of 
the  system  design  of  the  compound  existing  multi-variable  function 
curve  such  as  this  type  lies  in  what  future  single  variable  signals 
(x.j,  x^,...,  xn)  of  its  own  sensor  to  synthetically  merge  into  con¬ 
trol  signal  f(x^,  X2>...»  xn)  that  are  in  accord  with  curve  requirements, 
and  that  these  function  curves  are  obtained  in  model  design  phases 
through  wind  tunnel  tests  and  computations.  In  order  to  satisfy  the 
requirements  through  test  flights  seeking  the  most  favorable  curve,  the 
design  of  the  system  should  be  able  to  vary  the  predetermined  curves 
within  a  fairly  large  range. 

Using  differential  treatment  aircraft  to  realize  control  of  the 
multiple  element  function  is  feasible,  but  establishing  contact  with 
the  aircraft  load  and  adjusting  the  curve  is  difficult,  therefore  we 
can  only  seek  the  solution  computation  device  for  an  aircraft  electronic 
simulated  test.  The  methods  of  the  cone  model  cam  of  the  adopted  mach¬ 
inery,  the  electronic  diode  function  generator,  the  electronic  profile 
model,  etc.  to  realize  the  two-variable  function  are  comparatively 
complicated:  "Using  an  electronic  method  to  resolve  binary,  function 
still  up  to  now  has  had  almost  completely  no  development"  *  .  "The  re¬ 
search  work  of  reproducing  two  or  even  more  variable  function  trans¬ 
formation  devices  still  has  not  even  remotely  reached  the  extent  of 
perfection,"  "Three  or  even  more  variable  functions  using  the  above 
method  simply  cannot  be  reproduced" .  It  is  thus  evident  that  these 
methods  used  to  control  systems  are  difficult,  so  far  as  to  adjust 
transforming  two-variable  functions  then  is  even  more  difficult  to  im¬ 
agine.  This  text  proposes  a  method  for  simplified  engineered  realization 
for  some  separable  multi-variable  functions  to  become  single-variable 
functions,  and  designs  a  signal  synthe'ic  bridge.  This  not  only  resolves 
the  realization  and  adjustment  of  a  two-variable  function  bypass  door, 
but  also  possesses  reference  value  for  the  multi-variable  function 


II.  THE  REALIZATION  OF  THE  BYPASS  DOOR  TWO-VARIABLE  FUNCTION  CURVES 

Certain  high  speed  aircraft  bypass  door  system  requirements  are 
according  to  the  two  parameters  of  total  temperature  and  M  number  to 
carry  out  control,  its  open  angle  is  within  the  range  of  Cf~— 10”  with 
the  rise  and  increase  of  total  temperature,  and  at  the  same  time  with 
the  increase  and  reduction  of  M  number.  Its  control  curves  are  shown 
in  Fig.  1.  In  the  figure  L  is  the  effective  open  degree  of  the  bypass 
door;  S  is  the  rectilinear  displacement  of  the  control  motion  .tube. 
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FIG.  1:  CONTROL  CURVES  OF 
THE  BYPASS  DOOR 


1 .  RESOLUTION  OF  THE  BYPASS  DOOR  TWO-VARIABLE  FUNCTION  CONTROL 
CURVES 

In  Fig.  1  the  mathematical  expression  of  the  two-variable  curve  is: 

(  0  T<T , 

Sm  <  r,<r<r;  n-r.+rj-r, 

(  19.1*  t>ti  d) 


T,+150(M  — Mi) 
I  T. 


M<M. 

M>M, 


This  way,  a  curve  S(T,  M)  in  Fig.  1  can  resolve  as  the  two  nonlinear 
curves  that  Fig.  2  shows,  but  they  are  not  independent,  because  Tq  is 
the  function  of  the  M  number.  The  aircraft  load  profile  formula  atmo¬ 
sphere  data  computer  (central  apparatus)  only  can  export  the  single¬ 
variable  function  signal  of  the  M  number  or  total  temperature.  It 
cannot  give  out  extensive  function  f(  T-Tq(M)),  therefore  it  is  un¬ 
necessary  with  S(  T,  M)  to  separate  into  two  single  functions  S^(T) 
and  S2(M)  that  are  not  interrelated. 


FIG.  2:  DIVISION  OF  THE  CURVES 
IN  FIG.  1 

2.  CONDITIONS  FOR  SEPARATION  OF  MULTI-VARIABLE  FUNCTIONS 


By  unfolding  the  equation  of  the  binary  function  we  can  know  that 
so  long  as  their  intersecting  incline  derivative  sum  is  zero,  if: 


*5(7%  M) 


0 


(3) 


then  two  non-interrelated  single-variable  functions  that  are  separable 
are  engineered.  If  tS(T,  M)/?T  is  not  related  with  the  M  number, 
gS(T,  M)/aM  is  not  interrelated  with  T,  then  formula  (3)  is  bound  to 
obtain  satisfaction.  These  kinds  of  curves  inevitably  are  parallel  with 
one  another,  if  a  certain  basic  curve  follows  the  parallel  displacement 
of  a  certain  coordinate  axis,  the  distance  between  the  curves  can  be 
unequal,  but  they  must  be  parallel.  This  is  the  fully  required  condition 
for  a  two-variable  function  to  be  able  to  carry  out  separation.  It  can 
extend  to  multi-variable  functions: 


Arbitrary  multi-variable  function  f(x^,  x2,...,  xn),  if  its 
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second  intersecting  inclined  derivative  completely  is  zero, 

*/(■*!>  ”’j  *£>  _  0  ,  i  ,  j  —  1 ,  2  ,  * 

<  *  i  ( 4 ) 

thus 

/(*„  -♦-/*(■*»>  — •-/•<•*•> 

This  then  is  the  theoretical  basis  of  the  engineered  realization  of  a 
multi-variable  function  control  system. 


As  to  the  binary  function  is  Fig.  1,  within  the  scope  of  boundary 
condition  limitations  is  the  linearized  straight  line  of  the  distance 
equally  between  the  curves,  it  can  separate  as  two  linear  functions, 
and  is  a  particular  condition.  If  it  is  a  curve  where  the  distance  be¬ 
tween  the  curves  is  not  equal  but  is  mutually  parallel,  then  it  can 
separate  as  two  single-variable  nonlinear  functions  that  are  not  inter¬ 
related.  Taking  T=T^  and  M=M^  as  primary  points,  with  the  control 
curves  in  Fig.  1  the  region  between  Ty^T^T^,  and  develops,  if: 


Boundary 


conditions 


S(T,  M)-S,(r)+S.(M) 

(6) 

1 

r  o 

rcr, 

S,(T)-  - 

!  o.3688(r-r.) 

r,<T<Ti 

(7) 

i 

C  41.31 

T>T'% 

i 

r  o 

M<M, 

-56.32CM-M,) 

(8) 

(  -22.13 

M>M. 

0  <  S  (  r ,  M ) r  )  +  S,( M X 19 . 18 

(9) 

By  the  control  motion  tube  and  flight  enclosed  line  there  are  limitations. 


3.  THE  ENGINEERED  REALIZATION  OF  THE  BYPASS  DOOR  TWO-VARIABLE 
CONTROL  CURVE 


Using  the  electric  profile  model  to  reproduce  formula  (7)  and  (8) 
(See  Fig.  3),  with  resistance  characteristics  as: 
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FIG.  3:  SEPARATION  OF  THE  CONTROL 
CURVES  OF  THE  BYPASS  DOOR  AND  THE 
RESISTANCE  CHARACTERISTICS  OF  R(T) 

AND  ?.(M) 

j  R(T)~r.S,(T ) 

1  /?(M)-r.Sl(M)  +  22.13r.  (10) 

The  potentiometer  separately1  is  located  in  the  total  temperature 
sensor  and  the  M  number  computer,  when  the  control  apparatus  puts  out 
resistance  transformation  rc  ohms,  the  control  motion  tube  shifts  1mm. 
Substituting  formula  (6)  with  formula  (10)  we  obtain: 

S(Tt  M)“~-C /tf ( D+  —22.13  (11) 

In  order  to  solve  the  curve  of  formula  (11)  realized  in  Fig.  1,  we 
designed  the  bridge  circuit  that  Fig.  4  shows,  with  the  total  temper¬ 
ature  exhaust  potentiometer  R^,  and  M  number  exhaust  potentiometer  R^ 
strung  together  on  bridge  arm  ob,  its  linking  equation  is:  If  total 
temperature  rises  then  R(T)  increases;  if  M  number  increases  then  R(M) 
reduces.  By  formulas  (10)  and  (6)  we  can  obtain  ob  supporting  circuit 
total  resistance: 


The  left  and  right  door  linear  return  potentiometer  exhaust  value  R(s) 
and  motion  tube  movement  distance  S  are  in  a  direct  ratio,  if: 

S(«)-r,S  ^3) 
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FIG.  4:  SYNTHETI'C  BRIDGE  OF  THE  CONTROL 

SYSTEM  FOR  THE  BYPASS  DOOR 

KEY:  (a)  Amplifier;  (b)  Servo  valve;  (c) 

Motion  tube;  (d)  Right  bypass  door;  (e) 

Amplifier;  (f)  Servo  valve;  (g)  Left  by¬ 
pass  door. 

Where  rQ  is  the  resistance  on  the  return  potentiometer  unit  length. 
Taking  R^=R^=22. 1 3tq ,  then  cb  or  db  support  loop  total  resistance: 

22.13r,  (14) 

causes  the  system  to  possess  a  fairly  high  degree  of  exhaust  sensitivity 
under  standard  atmosphere  temperature  within  the  range  of  M^-SM<M^  (See 
Fig.  1),  taking: 

i?i-3Qr„  ^  (1 5) 

Rq  is  the  bridge  piezoelectric  resistance,  in  order  to  regulate  the 
feedback  system  of  the  return  exhaust  potentiometer,  and  through  system 
coordinate  computations,  take  Rq=500  Q  . 

When  the  total  temperature  T  or  M  number  change,  the  unequal  piezo¬ 
electricity  between  oc(od)  after  being  corrected  through  amplification 
controls  the  servo  valve,  drives  the  motion  tube,  controls  the  bypass 
door,  and  at  the  same  time  changes  the  return  resistance.  Due  to  the 
small  valve  ring,  the  frequency  band  of  the  hydraulic  servo  system  is 
5^10  cycle/second,  and  flight  M  number  and  total  temperature  are  long 
cycle,  slow  transformation  signs.  Therefore,  the  left  and  right  bypass 
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doors  throughout  conduct  regulation  of  two-variable  function  curve 
according  to  what  is  shown  in  Fig.  1.  Now  taking  the  left  door  as  an 
example  we  will  prove  whether  or  not  it  reproduces  curve  requirements. 
Because  the  bridge  throughout  is  situated  in  balanced  conditions,  there¬ 
fore  we  have: 

(i6) 

With  R^  ,  R^,  and  R^  substituting  in  the  above  formula  and  by  formula 
(13)  we  can  obtain: 

S~~-CR(T)+  J?(M))-22.13  (17) 

' • 

The  right  side  of  the  equal  sign  of  the  above  formula  and  formula  (11) 
are  completely  identical,  this  shows  the  distance  S  of  the  movement  of 
the  motion  tube  is  in  complete  accord  with  the  binary  function  curve 
control  S(T,  M)  that  Fig.  1  determined. 

III.  ADJUSTMENT  OF  THE  TWO-VARIABLE  FUNCTION  CURVE  OF  THE  BYPASS  DOOR 

In  order  to  satisfy  the  requirements  through  a  test  flight  seeking 
the  best  curve,  we  can  manufacture  a  certain  number  of  series  of  in¬ 
stallations  central  apparatus  components  that  have  different  control 
curves  to  provide  for  flight  test- selection.  However,  the  central 
apparatus  puts  out  ten  kinds  of  control  signals  towards  each  system 
component  of  the  whole  aircraft.  Therefore,  in  reality  it  is  very 
difficult  to  satisfy  this  requirement  of  the  system  in  this  text,  and 
having  to  transform  the  output  component  inside  the  central  apparatus, 
then  the  end  product  is  costly,  and  adjustment  is  not  convenient.  The 
adjustable  resistance  R^Rj.  that  is  set  up  on  each  arm  of  the  synthetic 
bridge  of  this  aircraft  bypass  door  system  not  only  can  compensate  for 
the  null  deviation  of  the  left  and  right  return  transport  potentiometer, 
what  is  even  more  important  is  through  them  t  ie  binary  function  curve  \ 
curves  of  the  central  apparatus  can  be  conveniently  turned  and  shifted. 

1.  The  Influence  of  the  Change  in  R1  on  the  Two-Variable  Function 
Control  Curves 

Assume  the  bridge  early  matched  condition  is  R10=30rc,  R20=R40=^0ro ' 
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and  R^q=R^q=22. 1 3rQ .  When  R^  by  early  condition  R^g  transforms  - 
if  R^=R^g+AR^,  the  theoretical  curve  of  the  bypass  door  motion  tube  is: 


S 


1  r  r(t 2 

r.  L  1  +  )*, 


R(  M)  1 
1  +  A5,  J 


22.13 


(18) 


Contrasting  formulas  (17)  and  (18)  we  can  see  that  the  change  in  R^ 
corresponds  to  the  output  characteristics  that  change  the  total  pres¬ 
sure  sensor  and  M  number  computer,  the  equivalent  output  characterist¬ 
ic  equation  of  the  central  apparatus  after  the  change  in  is: 


R'(T ) 


R(T) 

1  +  Ar, 


,  R'  (M) 


RC  M) 
l  +Ar, 


(19) 


where  is  the  corresponding  value  of  the  change  in  R^  ,  if  AR-j  /R-)  q * 

The  influence  of  the  change  in  R^  is  shown  in  Fig.  5.  In  the  figure 
S  only  expresses  the  theoretical  required  value,  in  actuality  the  dis¬ 
tance  receives  the  nonlinear  limitations  of  the  motion  tube,  all  the 
negative  values  can  only  correspond  to  the  bypass  door  closed  position. 
(This)  can  prove  that  the  change  in  R^  fully  changes  the  beginning  tone 
pitch  Sq  and  the  end  tone  pitch  of  the  motion  tube  (if  a  prescribed 
minimum  and  prescribed  maximum  of  the  curve),  but  its  effective  distance 
with  the  product  of  R^  is  a  constant: 


R,'(.Si—St)  —  const 


(20) 


This  way  then  we  can  change  the  controlled  range  of  the  bypass  door 
through  adjusting  R^ ,  as  to  the  shifting  of  the  control  position  caused 
by  the  change  in  R1 ,  we  can  higher  and  lower  translate  the  curve  through 
R^  to  compensate. 


2.  The  Influence  of  the  Change  in  R„  on  the  Two-Variable 
Function  Control  Curve 


The  change  in  R2  also  can  equivalently  convert  the  change  for  the 
central  apparatus  output  characteristics.  Assume  R2  has  an  increment, 
for  example  ^2^20^  "t^2^  *  Then  we  have: 

j*.  -L-CO  +A£.mr)  +  (1  +A£_.mM)3-22.13  (20) 
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FIG.  5:  EFFECT  OF  CHANGE 
IN  R  ON  THE  CONTROL 

CURVES  OF  THE  BYPASS  DOOR 

In  that  case  the  equivalent  output  characteristic  of  the  central 
apparatus  is: 

+M,)R(T),  /?'(M)-(1  +A5,)/?(M)  (21) 

The  influence  of  the  change  in  Rg  is  shown  in  Fig.  6.  Comparing  it  to 
Fig.  5,  it  is  evident  that  the  increased  Rg  and  the  reduced  R 1  have 
equal  effects,  provided  1 +ARg=1 / ( 1 +AR^ ) ,  then  the  curves  they  correspond 
to  are  completely  identical. 


FIG.  6:  EFFECT  OF  CHANGE 
IN  Rp  ON  THE  CONTROL  CURVES 
OF  THE  BYPASS  DOOR 


3.  The  Influence  of  the  Change  in  R,  on  the  Two  Function  Control 
Curves  J 


By  Fig.  5  and  Fig.  6  we  can  see,  the  change  in  R^  or  Rg  can  change 
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the  effective  control  range  of  the  slope  and  bypass  door  of  the  binary 
function  curves,  but  with  bringing  the  bypass  door  up  to  the  shifting 
of  the  control  position,  in  order  to  make  the  control  angle  degree 
become  zero,  it  is  necessary  to  take  the  curves  after  adjustment  through 
the  changes  in  R^  and  and  move  entirely  along  the  vertical  axis 
from  top  to  bottom.  This  function  depends  on  controlling  R^  to  be 
realized. 

When  R^R^q+AR^,  by  the  bridge  balanced  condition  we  can  obtain: 

S-~CX(T)+J?(M)1-22.1.< — (22) 


DIG .  7:  EFFECT  OF  CHANGE 
in  r3  ON  THE  CONTROL  CURVES 

OF  THE  BYPASS  DOOR 

The  influence  of  the  change  in  R^  is  expressed  in  Fig.  7.  By  Fig.  7  and 
comparison  of  formulas  (22)  and  (17)  we  can  see  that  the  change  in  R^ 
does  not  change  the  condition  of  the  curve  and  makes  the  whole  curve 
up  and  down  shift  parallel.  When  R^  increases,  the  whole  curve  shifts 
down,  its  translational  quantity  is  equal  to  the  increment  of  R^  divided 
by  the  unit  length  resistance  value  of  the  return  transport  potentio¬ 
meter  . 

4.  Adjustment  Method  for  Control  Curves 

Simultaneously  changing  adjustable  resistance  R1 ,  R2»  and  R^ 
makes : 


* 
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(23) 


(•ffi-tfj.+AJ?,-.  f\.C  1  +  AS.) 
j*.-/?„  +  A *,-/?„(  i  +Ar,) 

( •/?!"/?„+ A/?, 

By  the  balanced  condition  of  the  bridge  we  can  derive  the  movement 
curve  of  the  motion  tube: 


(24) 


Substituting  formula  (24)  with  formula  (10)  we  obtain: 

S=kS(T,  M)+Sq  (25) 

Where : 


r.  ..  1  +Ab. 

1  +  AS, 

5. ’■22.13(4  -  1  )- 


Ai?. 


(26) 


Formula  (25)  shows,  due  to  the  changes  in  R1 ,  R2»  and  R3>  causes  the 
binary  function  curve  of  Fig.  1  to  amplify  the  back  section  of  k  com¬ 
pletely  along  the  vertical  axis  direction  translating  up  and  down  SQmm. 
With  the  boundary  conditions  formula  (9)  of  S(T,  M)  substitute  formula 
(25),  we  obtain  the  effective  control  range  of  the  motion  tube: 

AS-S.-S.-19.18  4  (27) 

As  to  the  arbitrarily  determined  bypass  door  control  range,  we  deter¬ 
mine  the  k  value  by  formula  (27).  As  to  the  already  known  k  value, 
based  on  formula  (26),  the  selection  of  and  R2  can  have  many  kinds, 
to  make  R^  and  R^  both  possess  comparatively  small  amounts  of  changes, 
assume : 

(AS, -AS  (AS,- AS 

A<1.  1  _  *>1.  i 

Iar,-  -AS|  Ia/?,--AS  (28) 

The  selection  of  An  seen  in  Fig.  8,  when  k<Tl  ,  we  can  f'voose  l/(l+AR,)» 

1  -AS,,  (  U— AS«)/(  1  +  ASi)i  when  k>1  ,  we  can  choose  i  +ASi,  1/(1— ..'Si), 

(  1 +Ar.)/(  1 -AS.)  .  Generally  we  select  AR  =  0.1-0.5.  This  way  it  is 
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both  convenient  to  adjust,  and  ensures  the  degree  of  sensitivity  of 
the  output  of  the  bridge. 


5.  Practicle  Examples  of  Adjustment 

Assuming  the  central  apparatus  determined  curve  as  shown  in  Fig.  1, 
the  required  design  adjusted  bridge  arm  resistance  makes  the  bypass 
door  automatically  controlled  according  to  the  curve  in  Fig.  9  within 
the  range  of  0"^-5°  according  to  the  two  parameters  of  total  temper¬ 
ature  and  M  number. 

The  5°  corner  of  the  bypass  door  through  the  rocker  arm  corresponds 
to  the  rectilinear  displacement  =10. 92mm  of  the  motion  tube,  requiring 
Sq=0,  therefore: 

*,S  =  Si  —  S,~  10.92mm  (29) 

By  formula  (27)  we  obtain: 

19^8  *°-5693  (3C) 


■s<mm> 


io4 


o 


light  enclosed  line 


tahdard  temperature 
line 


FIG.  8:  EFFECT  OF  CHANGE 
IN  ON  k 


FIG.  9:  CONTROL  CURVES  OF  THE 
BYPASS  DOOR  FROM  0C to  5° 


Refering  to  Fig.  8,  taking  the  form  of  k -( 1  —  Ar*)/(  *  +  ^Ri)  (other  forms 
will  make  AR  enlarge),  and  substituting  the  already  known  k  value  we 
obtain: 

As  “0.2(45  (31) 

By  formula  (23)  we  obtain: 

£,  =  (<?„(  1  +  \R)  =  38.235r.  (32) 

R.  =  R.t(  1  -  \'r)=- 21.765  r,  (33) 

Making  Sq=0  in  formula  (26),  we  obtain  AR^,  and  already  known  4R^q=22. 
•  13r  ,  therefore  we  obtain: 

22. 13Ar,  =  12.06  r,  (34) 


the  right  door  parameters:  ^5=^30 

IV.  CONCLUSION 


The  control  target  of  many  flight  control  systems  except  carry¬ 
ing  out  control  according  to  a  direct  variable,  often  still  must  be 
affected  by  one  or  two  other  parameters  to  conduct  correction  and 
compensation.  This  type  of  designed  channel  of  open  formula  control 
system  of  multiple  intake-single  exhaust  takes  the  multi-parameter 
function  signal  synthesis  of  its  own  variable  sensor  as  a  control 
signal,  it  later  is  the  design  and  synthesis  of  the  single  intake- 
single  output  system  of  constant  curves.  The  method  of  engineered 
realization  of  separation  of  the  two-variable  function  control  curve  of 
the  bypass  door  that  this  text  expounds  on,  through  practicle  proof  is 
feasible.  This  method  can  extend  to  the  multi  element  function  research. 
The  engineered  multi-variable  function,  many  permitted  similar  simp¬ 
lifications  are  separable  multi-variable  functions.  Their  mathematical 
weakness  is  their  two  intersecting  partial  derivatives  completely  are 
null,  the  properties  in  the  figure  are  parallel  shifts  of  one  one-dim¬ 
ension  control  function  in  n  dimensional  space.  This  type  of  function 
can  become  a  single-variable  function  sought  to  be  realized. 

Through  the  method  of  adjusting  the  bridge  arm  resistance  to  change 
the  binary  function  control  curves,  installation  control  test  proof  is 
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■S 


simply  and  conveniently  feasible.  It  not  only  makes  the  bypass  door 
curve  along  the  vertical  axis  direction  adjustment,  but  also  can,  within 
a  specified  range,  make  the  whole  curve  along  the  model  axis  shift  to¬ 
wards  the  right.  This  is  through  adjusting  through  and  making 
S^O  in  formula  (26),  and  the  realistic  movement  of  the  motion  tube 
can  only  realize  the  S>0  section  to  be  completed.  By  Fig.  1  we  can 
know,  along  the  model  axis  the  largest  shiftable  range  is AT=T^-T^. 

Aside  from  this  we  still  must  point  out,  in  on-site  aircraft  tests,  we 
only  must  grasp  the  influence  of  the  changes  in  R^  through  R^  on  the 
control  curves  and  use  the  collective  test  method  of  controlling  while 
measuring.  We  do  not  use  the  theoretic  computation  if  we  can  make  the 
bypass  door  curve  based  on  necessary  conducted  adjustment. 

The  method  for  an  engineered  realization  of  the  separable  multi- 
variable  function  that  this  text  introduces,  and  the  methods  of  the 
synthetic  bridge  of  Fig.  4  and  the  bridge  arm  resistance  going  through 
adjustment  to  change  the  curve,  possess  a  specified  reference  value 
to  the  research  of  analogue  devices  of  the  compound  existing  three 
above  variable  functions  and  the  engineered  design  of  a  multiple  in¬ 
put-single  exhaust  system. 
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STRAIN  FATIGUE  PROPERTIES  AND  FATIGUE  DISLOCATION  STRUCTURES  OF 
2024  ALUMINUM  ALLOY 


Li  Qingsheng,  Tsai  Chikung,  Chen  Xianxi,  Pan  Tianxi 
Central  Iron  and  Steel  Research  Institute 
Zhou  Zhaoxiang 
Hongan  Corporation 

Abstract : 


This  paper  mainly  deals  with  the  relation  between  the  strain  fatigue  pro¬ 
perties  and  the  fatigue  dislocation  structures  in  2024  aluminum  alloy.  Some 
interesting  features  of  the  fatigue  dislocation  structure  in  this  commercial  alloy 
are  observed. 

The  analysis  of  deformation  work  density  for  smooth  fatigue  speci¬ 

mens  shows  that  the  strain  fatigue  properties  of  specimens  aged  at  240”C  are 
much  better  than  naturally  aged  specimens  at  high  strain  level. 

Systematic  study  of  fatigue  dislocation  structures  in  this  alloy  under  the 
transmission  electron  microscope  reveals  that  wavy  slip  mode  at  high  strain  . 
level  and  plain  slip  mode  at  low  strain  level  are  characteristics  of  the  dislo¬ 
cation  slip  in  the  naturally  aged  specimens,  in  contrast  cross  slip  is  much  easier 
ij  the  240*C  aged  specimens  at  different  strain  levels.  By  foil  thickness  ana¬ 
lysis  the  existence  of  {112}  slip  plane  in  240*C  aged  specimens  was  confirmed. 

It  is  also  pointed  out  that  the  helical  dislocations  in  the  alloy  are  formed 
through  cutting  of  moving  dislocations  with  dislocation  loops. 

The  conclusion  has  been  drawn  that  the  features  of  dislocation  slip  serve 
as  an  essential  factor  which  affects  the  interaction  of  dislocations  with  other 
microstructures  and  determines  the  fatigue  properties. 

I.  MATERIALS  AND  ACTUAL  TEST  METHODS 

We  used  two  kinds  of  aging  system  smooth  cylinder  specimens  to 
carry  out  strain  control  fatigue  tests,  and  separately  researched  the 
dislocation  structure  of  forward  added  load  and  rear  fatigue  cracks. 
The  material  selected  is  standard  component  LY-12  aluminum  alloy. 

The  two  kinds  of  aging  systems  are: 

Naturally  aged:  500cC,  two  hr  water  tempered,  aged  7  days  or  more. 
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Artificially  aged:  500° C,  2  hr  water  tempered,  aged  230* C  4  hr. 


The  one-way  tensile  yielding  strength  of  the  two  kinds  of  aging 
system  specimens,  the  tensile  strength,  and  the  cycle  yielding  strength 
are  seen  in  Table  1 . 

TABLE  1 .  MECHANICAL  PROPERTIES  OF  NATURALLY  AGED  AND  240  C*  AGED  SPECIMENS 


Yielding 

Strength 

a. 

Tensile 

Strength 

<j* 

Elasticity 

E 

Cycle  Yielding 

Strength 

Naturally  Aged 

42.5 

59 

7500 

47.5 

Artificially  Aged 

32 

42 

7600 

28 

(units  are  all  kg/rnm2-) 


The  specimen  is  a  6  13  round  stick  sample,  the  surface  has  a  smooth 
finish  V8  or  higher.  Tne  specimen  conducted  strain  control  fatigue 
tests  on  an  Instron  1255  test  computer.  After  the  cycle  crack,  at  the 
source  spot  of  the  crack  opening,  vertical  to  the  crack  surface  using 
a  line  corresponding  to  the  computer  we  take  a  0.3mm  piece  of  foil.  For 
teh  specimen  of  no  added  load,  we  take  a  piece  of  foil  of  the  same  thick¬ 
ness  from  the  specimen  central  section.  After  using  up  to  600#  sand 
paper  and  polishing  up  to  50  ju  with  a  foil  piece,  we  use  a  solution  of 
2:1  methanol:  nitric  acid,  and  a  dual  injection  electrolysis  buff.  The 
foil  is  observed  in  a  J3EM-200  electron  microscope. 

II.  RESULTS  AND  ANALYSIS  OF  ACTUAL  TESTS 

1 .  The  strain  fatigue  actual  test  results  and  their  method  of  de¬ 
formation  work  degree  are  described.  Test  results  of  strain  fatigue 
with  the  specimens  of  the  two  kinds  of  aging  treatments  and  the  use  of 
molded  strain  amplitude  -i— are  described,  as  shown  in  Fig.  1.  At 
high  strain  level,  the  relationship  between — phe,  the  artificially  aged 
specimen  and  fatigue  life  Nf  can  be  expressed  as: 

03  (1) 
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FIG.  1:  RELATION  BETWEEN 
AND  Nf 

'Aged  at  240°  G 
+Aged  naturally 

The  naturally  aged  specimen  can  be  expressed  as: 

-|-A*,-0.003  (2) 

3  L 

In  the  range  of  low  strain  (life  reaches  lO-^IO4’),  the  test  data 
all  clearly  diverges  the  Coff in-Manson  model  relationship  such  as 
formulas  (1)  and  (2)  express. 

Literature  flj  and  (2)  separately,  from  theory  and  actual  tests, 
pose  and  prove  within  the  range  of  extensive  life  (10^-10^  times), 
the  use  of  deformation  work  degree  to  replace  the  molded  strain  ampli¬ 
tude,  and  regard  as  strain  fatigue  damage  to  describe  the  rationalization 
of  the  parameters. 

We  use  the  deformation  work  degree  to  sort  out  the  data  of  Fig.  1. 
The  results  are  as  shown  in  Fig.  2. 

The  relation  between  the  artificially  aged  specimen  deformation 
work  degree Jadtu  and  Nf  can  be  expressed  as: 

|od«-3.22N;*-,M  (3) 
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FIG.  2:  RELATION  BETWEEN  °d. 
AND  Nf  J 

•Aged  at  240*  C 
+Aged  naturally 


In  accordance  with  the  equation  for  the  naturally  aged  specimen: 


"J"  odt  =■  1 . 24  jV}0-**  (4) 

As  Fig.  2  shows,  the  slope  of  the  artificially  aged  specimen  is 
fairly  large,  following  the  reduction  of  deformation  work  degree,  the 
data  of  the  naturally  aged  specimen  is  gradually  approached.  Therefor 
artificially  aged  specimens  have  even  more  superior  strain  fatigue 
performance  when  enduring  comparatively  high  strain. 


FIG.  3:  FATIGUE  DISLOCATION  IN 
NATURALLY  AGED  SPECIMENS  AT  HIGH 
STRAIN  LEVEL  (Nf»15) 


2.  Strai  .  dislocation  structure  and  analysis.  Microscopic 
analysis  of  the  foil  transmission  reveals:  the  difference  of  the  fatigue 
dislocation  structure  of  the  two  kinds  of  aging  system  samples  at 
different  strain  levels  is  unusually  striking. 

The  naturally  aged  specimen:  In  the  range  of  high  strain,  dis¬ 
location  bends,  for  the  time  being  are  uniformly  distributed,  demon¬ 
strating  wavy  slip  characteristics.  Figure  3  shows  this  kind  of  char¬ 
acteristic,  in  the  lower  part  of  Fig.  3  is  the  inclusion  (Fe,-  Mn,  Si) 

Ai^  of  the  crack  opening  under  the  effect  of  dislocation.  In  this 
kind  of  condition,  microcosmic  mould  strain  distribution  is  uniform, 
the  fatigue  damage  distribution  of  the  inclusion  open  crack  formation 
is  also  uniform,  this  is  simnly  the  reason  for  the  macroscopic  quantity 

— j— and  fatigue  life  Nf  in  accord  with  the  Cof f in-Manson  relation. 

Following  the  reduction  of  strain  quantity,  the  degree  of  disloc¬ 
ation  reduces,  when  strain  quantity  is  less  than  the  strain  quantity 

3  L 

that  corresponds  to  10^10  number  of  cycles,  the  dislocation  config¬ 
uration  has  substantial  change. 

In  the  range  of  low  strain,  the  typical  dislocation  configuration 
is  as  shown  in  Fig.  4.  Figure  4-A  shows  the  gathering  and  remaining 
wavy  slip  band  of  inclusion  around  dislocation. 


FIG.  4:  FATIGUE  DISLOCATION, IN  NATURALLY  AGED  SPECIMENS  AT 
LOW  STRAIN  LEVEL  (Nf =1 . 2x1 05 ) .  A.  Electron  beam  direction  (112) 


Figure  4-B  shows  dislocation  of  dislocation  mesh  and  strategic  position 
area  of  inclusion  of  the  specimen  surface  form. 


Because  the  dislocation  plane  is  wavy  slip,  the  form  remains  a 
wavy  slip  band.  In  the  band  the  GP  region  repeatedly  sheared  and 
resolved  by  displacement,  this  kind  of  microcosmic  organization  de¬ 
teriorates,  and  promotes  the  fatigue  crack  veins  in  the  band  to  form 
and  expand  in  the  early  stage.  Because  dislocation  is  highly  con¬ 
centrated  in  the  wavy  slip  band,  microcosmic  strain  distribution  is 
not  uniform.  Moreover,  fatigue  damage  in  the  wavy  slip  band  that  the 
GP  region  helps  bring  about  also  is  local.  This  way,  using  macroscopic 
quantity  — 1  to  express  the  Cof f in-Manson  relation  of  fatigue 

damage  to  forecast  fatigue  life,  then  the  results  obtained  will  be  much 
higher  than  in  reality. 

Artificially  aged  specimen:  Figures  5  and  6  separately  express 
the  microcosmic  organization  and  fatigue  dislocation  structure  in  the 
range  of  high  and  low  strain.  As  Fig.  5  and  Fig.  6  show,  artificially 
aged  specimens  have  alike  foil  condition  S'  that  the  edge  base  of£j00^ 
plane  separates  out.  In  the  range  of  low  strain,  this  alike  S'  can 
have  advancement  towards  the  formation  of  fatigue  crack  veins*  creat¬ 
ing  divergence  to  the  Coff in-Manson  relation. 

As  Fig.  5  and  Fig.  6  show,  in  the  range  of  high  strain,  degree  of 
dislocation  is  high,  in  the  range  of  low  strain,  degree  of  dislocation 
is  low;  additionally,  on  the  line  of  dislocation  all  the  sediment  is 
separated  out. 


FIG.  5:  FATIGUE  DISLOCATION 
IN  ARTIFICIALLY  AGED  SPECIMEN 
AT  HIGH  STRAIN  LEVEL  (Nf=225) 


001  polar  circle 

FIG.  6:  POLE  FATIGUE  ANALYSIS 
OF  DISLOCATIONS  IN  (211)  PLANE 
KEY:  (a)  trace. 

With  artificially  aged  specimens,  dislocation  is  easy  to  interact 
with  wavy  slip.  Aside  from  dislocation  in  111  plane,  we  still  have 
dislocation  in  112  plane.  As  Fig.  6  shows,  three  comparatively  wide 
dislocations,  their  trace  lines  (the  connecting  lines  of  the  extreme 
points  of  dislocation)  are  110  direction,  pole  analysis  reveals,  they 
possibly  are  dislocations  in  plane  111.  The  upper  left  side  dislocation 
trace  lines  are  120  direction,  possibly  are  dislocations  in  211  plane, 
neighboring  dislocations,  and  should  correspond  to  the  same  foil  thick¬ 
ness.  Different  crystallography  planes  with  foil  plane  included  angle 
a,  degree  of  dislocation  s,  and  foil  thickness  t  should  s- tisfy: 

I  -*tfo  ( 5 ) 

Corresponding  to  neighboring  dislocation,  we  should  have: 


Plane  211  and  plane  001  (foil  plane)  included  angle  is  65.9°, 
plane  111  and  plane  001  included  angle  is  54.73°;  actual  measurement; 
of  s.j  is  6mm,  and  S£  is  4mm.  The  above  data  is  basically  in  accord 


with  formula  (6).  By  this  we  can  prove,  the  upper  left  dislocation 
of  Fig.  6  is  dislocation  in  the  211  plane. 

Literature  fi)  had  produced  the  existing  wavy  slip  112  plane  in 
pure  aluminum  60  C  tensile  specimen.  LY-12  aluminum  alloy  artificially 
aged  specimens  in  conditions  of  room  temperature  and  cyclic  load  have 
dislocations  existing  in  112  plane.  This  fact  points  out  the  possi¬ 
bility  of  wavy  slip  in  112  plane.  The  lower  left  dislocation  in  Fig. 

5,  possibly  is  the  condition  of  wavy  slip  intersecting  with  111  plane 
by  211  plane. 

With  artificially  aged  specimens,  dislocation  is  easy  to  inter¬ 
sect  with  wavy  slip.  This  way,  the  possibility  of  forming  a  remaining 
wavy  slip  band  and  a  strategic  area  of  inclusion  reduces.  This  is  the 
reason  for  the  favorable  strain  fatigue  performance  of  the  artificially- 
aged  specimen  compared  to  the  naturally  aged  specimen.  E.J.  Coyne's^' 
work  on  the  Al-Zn-Mg  alloy  also  produced  a  good  fatigue  performance 
of  aged  samples,  and  attributed  the  reason  to  the  dislocation  being 
easy  to  intesect  with  wavy  slip. 

Helical  dislocation:  Helical  dislocation  is  a  particular  kind  of 
dislocation  of  the  above  alloy.  Artificially  aged,  no  added  load  con¬ 
dition  then  has  a  certain  amount  of  helical  dislocation.  Figure  7  is 
the  appearance  of  the  dark  field  of  the  weak  bond  of  helical  dislocation 
under  this  condition.  A  shows  the  intersection  of  a  moving  dislocation 
and  a  hoop  dislocation.  The  lower  right  side  of  A  still  has  an  in¬ 
dependent  hoop  dislocation.  We  point  out  that  the  above  helical  dis¬ 
location  system  was  formed  through  cutting  of  moving  dislocations  by 
the  empty  position  collapsing  and  the  formed  hoop  dislocation  that  were 
produced  by  fire  tempering. 

In  artificially  aged  specimens  after  fatigue  cracks,  the  changes 
in  form  and  appearance  of  helical  dislocation  are  not  large. 

Naturally  aged  specimens  and  conditions  of  no  added  load  only  have 
very  small  amounts  of  helical  dislocation.  After  fatigue  cracks,  hel- 
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FIG.  7:  HELICAL  DISLOCATION 
IN  ARTIFICIALLY  AGED  SPECIMEN 
BEFORE  LOADING 
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FIG.  8:  HELICAL  DISLOCATION 
IN  NATURALLY  AGED  SPECIMEN 
AFTER  FATIGUE  CRACKS 


icai  dislocation  amount  increases.  Figure  8  shows  the  high  degree  of 
dislocation  concentration  that  helical  dislocation  produces  after  fat¬ 
igue  cracks.  In  artificially  aged  specimen  comparisons,  this  possibly 
is  related  to  the  intersection  of  wavy  slip  performance  with  dislocation. 

III.  CONCLUSION 

1 .  Using  the  degree  of  deformation  work  we  can  describe  very  well 
the  strain  fatigue  laws  of  different  time  conditions  with  the  LY-12 
aluminum  alloy.  In  the  range  of  high  strain,  artificially  aged  spec¬ 
imens  have  obvious  advantageous  strain  fatigue  performance. 

2.  In  heat  treatment  systems,  strain  quantities  all  change  fatigue 
dislocation  configurations  and  wavy  slip  characteristics. 

3.  Artificially  aged  specimens,  under  conditions  of  room  temper¬ 
ature  and  cycle  load,  have  dislocation  existing  in  112  plane,  and  can 
have  wavy  slip  in  112  plane. 

4.  Helical  dislocation  systems  in  LY-12  aluminum  alloy  are  formed 
by  cutting  of  moving  dislocations  with  dislocation  loops. 
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STUDY  ON  MICROSTRUCTURS  OF  CARBON-ALUMINUM  COMFOSITES 


Zheng  Lide,  Zhang  Guoding,  Wang  Wenrong,  Shi  Peiyuan 
Shanghai  Communications  University 

Abstract: 

The  microstrueture  of  carbon-aluminum  composites  is  close  related  with  their 
technological  process.  Therefore,  the  problems  in  the  process  can  be  discovered 
by  studying  the  features  of  the  microstructure,  and  available  information 
can  be  obtained  further  to  choose  the  process  parameters  correctly,  to  improve 
the  technological  process  and  to  enhance  the  properties  of  tht  fibers  and  compo¬ 
sites. 

The  microstructures  of  the  carbon  fibers,  the  carbon  fibers  coated  with  nickle, 
the  carbon  fibers  with  nickle  layer  after  vacuum  heat  treatment,  carbon-alumi¬ 
num  composite  wires  and  their  specimens  fabricated  by  hot  pressing  under  the 
different  technological  conditions  have  been  taken  with  the  metal lographic  and 
sweep  electron  microscopes  and  are  presented  in  this  paper. 

• 

Carbon-aluminum  composite  material  is  a  new  type  of  composition 
of  material  that  has  had  rapid  development  in  the  last  twenty  years. 

It  has  a  higher  strength,  higher  mould  capacity,  high  temperature 
performance  that  is  good,  and  excellent  properties  of  dimensional 
stability.  This  material  has  very  great  potential  in  applications 
to  space  navigation  and  aeronautics  industries.  The  properties  of 
carbon-aluminum  composites  are  closely  related  to  the  technological 
processes  of  basic  type,  carbon  fiber  surface  coating  and  metal  liquid 
state  soaking  and  hot  press  solidification,  and  especially  are  directly 
related  to  the  microstructures  and  defects.  Therefore,  to  research 
the  microstructures  and  defects  of  plated  carbon  fiber  coats  after 
coating  treatment  or  carbon  fibers  with  basic  application,  and  carbon- 
aluminum  composites  after  soaking  and  hot  press  solidification  can  be 
to  research  a  definite,  rational  technological  process  and  parameters 
and  their  influence  on  the  microstructure  and  properties  of  carbon- 
aluminum  composites;  the  combining  and  interaction  between  fiber  and 

basic  (or  plated  coating)  provides  an  important  reference  and  foundation. 


150 


This  text  uses  the  method  of  an  optical  microscope  and  sweep 
electron  microscope  to  conduct  analysis  and  research  on  the  micro¬ 
structure  and  fracture  of  nickel-plated  carbon  fibers,  the  nickel- 
plated  carbon  fibers  after  vacuum  heat  treatment,  carbon-aluminum 
composite  wires,  and  carbon-aluminum  composite  material. 

EXPERIMENTS 

The  analysis  of  microstructures  of  carbon-aluminum  composites 
primarily  researches  the  characteristics  of  the  microstructure  and  its 
transformation  of  carbon  fibers,  coatings  and  carbon-aluminum  composites 
in  each  technological  process  that  prepares  the  carbon-aluminum  com¬ 
posite.  In  the  text  we  adopt  metallographic  analysis  method  to  invest¬ 
igate  the  uniformity  of  the  carbon  fiber  coating,  the  interaction  of 
heat-treated  carbon  fibers  and  nickel  layer,  the  distribution  of  fibers 
in  hot-press  solidified  carbon-aluminum  composites,  the  binding  con¬ 
ditions  of  fiber  and  basic,  and  carbon-aluminum  composite  wire  composites. 
In  metallographic  analysis,  because  there  is  a  big  difference  in  the 
properties  of  various  metals  in  the  specimen,  therefore  the  process 
of  preparing  the  specimen  has  a  very  great  influence  on  whether  or  not 
the  microstructure  and  defects  are  correctly  shown.  The  prepared 
specimen  can  use  two  methods  of  preparation  which  are  using  6101  epoxy 
resin  and  ethylene  diamine  solidifier  according  to  appropriate  propor¬ 
tions  to  carry  out  casting,  or  using  bakelite  powder  to  hot-press. 

Additionally,  we  still  use  the  sweep  electron  microscope  to  observe 
and  analyze  the  interaction  of  fractures  and  carbon  fibers  with  the 
nickel  of  the  nickel-plated  carbon  fiber  carbon-aluminum  composite  wires, 
and  the  hot-press  solidified  carbon-aluminum  composite  plates. 

The  preparation  method  and  technological  process  conditions  of  the 

various  specimens  that  this  text  used  already  have  detailed  explanation 

i'l—3'1 

in  related  texts  '  ,  we  will  not  again  repeat  it  here. 

1 .  The  Microstructure  of  Nickel-Plated  Carbon  Fibers 

In  order  to  improve  the  binding  between  carbon  fibers  and  alum- 
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inum,  we  use  a  chemical  nickel-plating  method  on  the  carbon  fiber 
surface  to  plate  one  layer  of  nickel.  Tests  prove  that  it  has  an 

( 3) 

excellent  result  in  promoting  carbon  fiber  and  aluminum  composites  v  '  . 
The  condition  of  the  binding  of  the  coat  with  the  carbon  fiber,  and  the 
thickness  of  the  coat,  and  the  degree  of  uniformity  all  are  closely 
related  to  the  technological  parameters.  Under  conditions  of  approp¬ 
riate  technology  we  can  plate  a  thin  and  uniform  nickel  coat  micro- 
structure;  see  Fig.  1-a.  When  the  plating  time  increases,  the  thickness 
also  correspondingly  increases,  as  Fig.  1-b  shows.  The  cleanliness  of 
the  carbon  fiber  surface  is  the  necessary  condition  to  obtain  a  good 
coating,  in  incomplete  rinsing  the  fibers  and  the  coating  have  an 
obvious  disjointed  appearance,  even  complete  coats  will  not  stay  on, 
as  shown  in  Fig.  1-c  and  1-d.  The  surface  appearance  of  a  typical 
nickel-plated  carbon  fiber  is  shown  in  Fig.  1-e. 

2.  Microstructure  of  Nickel-Plated  Carbon  Fibers  After  Vacuum 
Heat  Treatment 

The  observance  of  microstructure  of  nickel-plated  carbon  fibers 
after  vacuum  heat  treatment  is  an  important  method  of  researching  the 
interaction  of  nickel  carbon.  Following  the  change  in  temperature  and 
time  of  vacuum  heat  treatment,  the  microstructure  of  nickel-plated 
carbon  fiber  has  a  series  of  evident  changes. 

After  heat  treatment  below  700® C,  the  microstructure  of  the  car¬ 
bon  fiber  has  no  evident  changes.  When  the  temperature  raises  to  750° C 
and  treatment  time  increases,  a  "surface"  structure  appears  on  the 
nickel-plated  carbon  fiber,  and  the  nickel  in  the  nickel-plate  layer 
assumes  a  needle  form  that  permeates  into  the  inner  part  of  the  carbon 
fiber,  as  Fig.  2-a  shows.  When  the  temperature  continually  rises,  the 
nickel  coating  thins  (even  is  eliminated)  the  nickel  begins  to  build 
up  in  the  inner  part  of  the  carbon  fiber.  The  microstructure  after 
treatment  at  800°  C  is  shown  in  Fig.  2-b. 

3.  The  Microstructure  of  Carbon-Aluminum  Composit  Wires 

Using  chemical  gas  deposit  method,  the  deposit  on  the  carbon  fiber 
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surface  is  one  titanium-boron  coating,  and  directly  soaking  aluminum 
is  an  effective  method  of  preparing  the  carbon-aluminum  composite  wires. 
But  the  control  of  the  gas  deposit  and  soaking  technological  parameters 
has  a  very  large  influence  on  the  infiltration  and  composites  of  car¬ 
bon  fibers  with  aluminum  liquid.  When  the  control  of  the  technological 
parameters  is  suitable,  the  carbon  fibers  and  aluminum  liquid  can  spon¬ 
taneously  infiltrate,  and  fiber  distribution  is  also  comparatively 
uniform,  as  Fig.  3-a  shows.  When  the  control  of  technological  parameters 
is  not  good,  then  it  cannot  have  a  good  soaking  interaction  with  the 
aluminum,  as  shown  in  Fig.  3-b.  In  the  figure  it  is  evident  that  in 
the  center  part  aluminum  liquid  has  no  immersion.  Figures  3-c  and  3-d 
are  photographs  of  fracture  of  carbon-aluminum  composite  wires.  From 
the  center  it  is  evident  that  when  carbon  fibers  produce  boundary  sur¬ 
face  reaction  with  aluminum,  fiber  and  aluminum  interaction  strengthens 
the  fiber  fracture  level,  and  there  is  no  appearance  of  pulling  out. 

When  the  interaction  of  fiber  and  basic  is  appropriate,  the  fibers  have 
a  slight  pulling  out,  as  Fig.  3-d  shows. 

4.  Microstructure  of  Carbon-Aluminum  Hot-Pressed  Solidified 
Composites 

Under  temperatures  that  are  close  tojthe  melting  point  of  aluminum, 
carbon  fiber  or  carbon-aluminum  composit  wire  that  will  go  through 
coating  treatment  can  obtain  carbon-aluminum  composite  plate  with  al¬ 
uminum  foil  hot-press  solidification.  Figure  4-a  is  a  photograph  of 
a  carbon-fiber  and  aluminum  hot-press  specimen  that  has  no  coating, 
demonstrating  that  carbon  fibers  and  aluminum  do  not  have  a  very  good 
interaction.  Figure  4-b  is  a  photograph  of  a  carbon  fiber  and  hot-press 
specimen  that  has  a  coating,  demonstrating  that  interaction  of  carbon 
fibers  and  aluminum  basic  is  very  good.  Figure  4-c  shows  the  separation 
condition  of  carbon  fibers  in  the  carbon-aluminum  composit  wire  and 
aluminum  hot-press  specimen. 

CONCLUDING  REMARKS 

The  analysis  on  nickel  plated  carbon  fibers,  vacuum  heat  treated 
nickel-plated  carbon  fibers,  carbon-aluminum  composite  wires,  and 
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hot-press  test  samples  make  clear,  their  microstructure  and  its 
technological  process  have  a  close  relationship,  and  therefore  we 
can  discover  the  existing  problem  in  the  technological  process  through 
analysis  and  research  of  the  microstructure,  provide  a  foundation  to 
further  correctly  choose  the  technological  parameters,  improve  the 
technological  process,  and  enhance  the  properties  of  fibers  and  com¬ 
posite  material. 

The  metallographic  method  and  sweep  electron  microscope  method  are 
effective  means  to  analyze  and  research  the  microstructures  of  carbon 
fibers,  nickel-plated  carbon  fibers,  vacuum  heat  treated  nickel-plated 
carbon  fibers,  carbon-aluminum  composite  wires,  and  hot  press  test 
samples,  but  it  is  necessary  to  interact  with  the  correct  manufacturing 
method  and  a  good  buff  grinding  technique,  then  we  can  really,  object- 
ivelt  show  the  form  of  microstructures. 

Regarding  the  effects  of  microstructures  of  carbon  fiber  coating 
and  carbon  fiber  and  aluminum  produced  coating  and  fiber,  the  effect 
of  basic,  etc.,  the  problems  are  pending  further  research.  At  the  same 
time  as  the  suitable  above  research  work,  we  still  must  explore  new 
manufacturing  methods  and  adopt  new  techniques. 


FIG.  1 -a:  MICROSTRUCTURE  OF  FIG.  1-b:  MICROSTRUCTURE  OF 
CARBON  FIBER  COATED  WITH  CARBON  FIBER  COATED  WITH 

NICKEL  (in  25  minutes)  NICKEL  (in  30  minutes) 
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FIG.  1-e:  TYPICAL  END 
VIEW  OF  CARBON  FIBER 
COATED  WITH  NICKEL 
■$■000  x 


FIG.  2-a :  200x 
MICROSTRUCTURE  OF  CARBON 
FIBER  COATED  WITH  NICKEL 
TREATED  AT  750*0  (in  1  hour) 


FIG.  2-bs  5000x  FIG.  3-a:  250x 

MICROSTRUCTURE  OF  CAR-  MICROSTRUCTURE  OF  CARBON-AL- 
BON  FIBER  COATED  WITH  UMINUM  COMPOSITE  WIRE 
NICKEL  TREATED  AT  800*C 
(in  1  hour) 


155 


m 


v* 


-\&i 


s* 


v 


FIG.  3-b!  250x 
MICROSTRUCTURE  OF  PARTIAL 
COMPOSITE  CARBON-ALUMINUM 
WIRE 
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FIG.  3-c:  20 OOx 

FRACTURE  OF  CARBON-ALUMINUM 

COMPOSITE  (Strong  Bonding) 
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FIG.  3-d:  1500x 
FRACTURE  OF  CARBON-ALUMINUM 
COMPOSITE  (Intermediate 
Bonding) 


FIG.  4-a:  1 500x 

MICROSTRUCTURE  OF  UNCOATED 
HEAT-PRESSES  CARBON  FIBERS 
ALUMINUM  COMPOSITE 
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FIG.  4-b:  500x  FIG.  4-c:  1 0Ox 

MICROSTRUCTURE  OF  COATED  MICROSTRUCTURE  OF  CARBON-ALUMINUM 
CARBON  FIBERS- ALUMINUM  COMPOSITE  (Hot-pressed  from  carbon 

OILS  COMPOSITE  (Hot-pressed)fiber-aluminum  composite  wires  and 

aluminum  foils) 
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FLOURISHING  DEVELOPMENTS  IN  ACADEMIC  ACTIVITIES  OF  CHINA  AERONAUTIC 
INSTITUTE 

In  the  China  Aeronautic  Institute  in  the  period  from  August  to 
October  1981,  academic  activities  were  perfectly  dynamic.  China 
Aeronautic  Institute  successively  with  China  Automation  Institute 
jointly  convened  "Control  System  Models  and  Real  Technological  Exchange 
Conference",  jointly  convened  "The  Third  National  Crack  Mechanics 
Symposium"  with  China  Mechanics  Institute  and  China  Machinery  Engineer¬ 
ing  Institute,  jointly  convened  "Titanium  and  Titanium  Alloy  Symposium" 
with  China  Metallurgy  Institute,  jointly  convened  "Engine  Combustion 
Symposium",  "Engine  Heat  Transmission  Symposium",  Impeller  Machinery 
Pneumatic  Heat  Mechanics  Symposium",  and  "Engineering  Heat  Mechanics 
Symposium"  with  China  Engineering  Heat  Physics  Institute.  In  addition, 
each  committee  of  related  disciplines  at  China  Aeronautics  Institute 
successively  convened  "Computer  Auxilary  Aircraft  Design-Auxilary 
Geometric  Design-Auxilary  Manufacturing  Symposium"  (simply  called 
CAD/CAGD/CAM  Symposium),  "Flight  Mechanics  and  Test  Flight  Symposium", 
"Applications  of  Modern  Control  Theory  and  Electronics  in  the  Area  of 
Aeronautics  and  Space  Navigation  Symposium",  "Non-Metallic  Material 
and  Technology  Symposium",  "Optics  Metallography  and  Transmission  El¬ 
ectron  Lens  Symposium",  "Aircraft  Technology  Symposium",  "Aeronautics 
Maintenence  Engineering  Specialty  Committee  Establishment  Conference 
and  Aeronautic  Maintenence  Engineering  Theory  Seminar",  "Special 
Problems  of  Composite  Material  Seminar",  and  others,  totalling  seventeen 
symposiums.  In  August  in  Beijing  they  also  jointly  held  "Composit 
Material  Short  Term  Training  Course"  with  Beijing  Aeronautics  Institute. 

The  above  academic . activities  all  are  academic  conferences  that 
belong  to  the  discipline  committees  or  specialty  study  groups.  The 
representatives  present  at  the  meeting  unanimously  considered  that  the 
convening  of  these  symposiums  was  perfectly  essential,  and  also  very 
timely.  For  some  meetings  it  was  the  first  time  held  in  the  last 
ten-plus  years.  Each  related  unit  and  numerous  experts,  professors,  and 

science  and  technology  personnel  enthusiastically  delivered  their  thes¬ 
es,  and  enthusiastically  supported  the  academic  activities.  The  research 
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papers  that  were  received  at  these  meetings  were  comparatively  many, 
the  content  was  considerably  rich,  the  related  specialty  learning 
scope  was  also  very  extensive.  From  the  theses  quality  and  level 
the  fields  were  clearly  at  a  higher  level  than  previously. 

These  academic  activities  possess  the  following  characteristics: 
Firstly,  they  are  alike  or  nearly  alike  academic  activities  that  are 
held  jointly  with  brother  institutes  and  have  specialty  content,  many 
famous  experts,  professors,  and  science  and  technology  personnel  from 
within  the  country  were  invited  who  participated  in  the  activities,  thus 
making  the  scope  of  the  symposiums  even  more  extensive.  At  the  same 
time  it  also  enhanced  the  quality  of  the  academic  activities  and  re¬ 
flected  the  present  level  of  science  of  our  country's  related  disciplines. 
Secondly,  these  academic  conferences,  aside  form  exchanging  the  research 
results  and  experience  of  recent  years  in  the  various  areas  of  scientif¬ 
ic  research,  production,  teaching,  etc.,  still  invited  related  experts 
to  the  conferences  to  make  the  summary  academic  report  of  technological 
developments  and  trends  inside  and  outside  the  country,  for  example: 
the  National  Computer  Department  Vice  Department  Head  Guo  Pingxing 
at  the  conference  of  application  of  modern  control  theory  and  electron¬ 
ics  in  areas  of  aeronautics  and  space  navigation  wrote  the  comprehensive 
report  "The  application  of  Mini  Computers  in  teh  Areas  of  Aeronautics 
and  Space  Flight";  Professor  Wen  Zhuanyuan  of  Beijing  Aeronautics  In- 
stitu^  at  the  "Control  Systems  Model  and  Real  Technology  Symposium" 
wrote  "A  Report  on  Going  to  America  to  Observe  and  Study  Ground  Flight 
Analogue  Computer";  Pro^es^or  Yang  Peng j i  of  Xi  Bei  Industrial  Univer¬ 
sity  at  the  "CAD/CAGD/CAM  Gonference"  made  the  academic  report  on  "7760 
CAD/CAM  System  Engineering",  making  the  representatives  at  the  conferen¬ 
ces  to  broader,  their  academic  outlook,  obtain  inspiration,  and  received 
a  warm  welcome  form  conference  representatives. 

In  the  above  academic  activities  the  committees  of  each  disciple 
at  the  same  time  held  work  conferences,  stressing  the  research  of  re¬ 
lated  academic  organizational  and  structural  problems.  They  established 

the  "Aeronautics  Maintenance  Engineering  Specialty  Committee",  the 
"Electron  Specialty  Committee",  and  the  "Automatic  Control  Specialty 
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Committee".  They  deliberated  on  forming  twenty-plus  disciple  study 
groups.  The  work  meetings  also  researched  the  arrangement  of  academic 
activities  for  each  disciple  from  now  on,  and  put  forward  plans  for 
activities  in  1982  and  preliminary  tentative  plans  for  activities  from 
1983  to  1985-  They  pointed  out  for  academic  activities  from  now  on: 

We  should  closely  combine  the  realistic  needs  of  our  country's  present 
four  modernizations  to  develop  academic  activities  of  some  special 
topics.  The  scale  of  the  meetings  can  be  somewhat  small,  the  discussion 
content  can  be  somewhat  deep,  and  interaction  of  scientific  research, 
production,  and  teaching  practices  can  become  even  closer. 


"JOURNAL  OF  AERONAUTICS"  BRIEF  RULES  FOR  SOLICITED  CONTRIBUTIONS 


"Journal  of  Aeronautics"  is  a  comprehensive  academic  publication 
sponsored  by  China  Aeronautic  Institute,  already  published  and  avail¬ 
able  inside  and  outside  the  country  in  1980.  This  publication  takes 
Marxism,  Leninism,  Mao  Zidong  thought  as  its  guide.  It  thoroughly 
carries  out  the  guiding  principle  of  "Let  a  hundred  flowers  bloom,  let 
a  hundred  schools  of  thought  contend",  develops  academic  exchange, 
and  promotes  the  modernization  of  space  navigation  science  and  technol¬ 
ogy* 


This  publication  primarily  publishes  scientific  research  results 
in  the  realm  of  aeronautics  and  space  navigation  science  and  technology. 
The  reader  target  is  primarily  specialists  engaged  in  scientific  research, 
teaching,  desgin,  production,  and  application  in  the  areas  of  aeronaut¬ 
ics  and  space  navigation  science  and  technology. 

I.  This  publication  publishes  the  following  various  contributions: 

1 .  New  academic  achievements  related  to  aeronautics  and  space 
navigation  science  and  technology,  including  aerodynamic  mechanics 
and  flight  mechanics,  aircraft  design  and  structure  strength,  aero¬ 
dynamic  devices,  electronic  technique  and  automatic  control,  academic 
discussion  of  produciton  in  the  area  of  technique  and  material,  and 
achievements  and  experiences  of  new  theory  and  new  technique  applications. 

2.  Technological  reading  notes,  reports,  and  summarized  discussions 
of  related  aeronautic  and  space  navigation  science  and  technology  that 
has  been  newly  created. 

3.  Explanation  and  discussion  of  books  and  literature  that  are 
related  to  aeronautics  and  space  navigation  science  and  technology. 

II.  Requirements  for  Contributions: 

1.  The  content  of  the  analysis,  computations,  test  data,  etc. 
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of  the  contributions  must  be  correct  and  without  error,  structure 
must  make  every  effort  to  be  close, j  arrangement  of  ideas  must  be  clear, 
grounds  of  argument  must  be  well-knit,  and  writing  style  must  be 
succinct  and  easy  to  understand.  Each  essay  should  take  5000  charact¬ 
ers  as  appropriate,  and  the  length  of  the  produced  dissertation  must 
not  exceed  10,000  characters  (including  illustrations  and  tables). 

2.  Contributions  use  16  evolution  ruled  paper.  A  fountain  pen 
must  be  used  for  writing.  Writing  must  be  neat  and  clear.  The  units  of 
technical  terms  and  computations  should  be  consistent  throughout. 

3.  Illustrations  shall  use  drawing  paper  according  to  the  engineer 
ed  drawing,  the  words  on  the  illustration  shall  be  written  with  a  pen¬ 
cil.  Please  position  illustrations  centered  on  the  manuscript,  and  use 
three  ruled  spaces  to  frame  the  drawing.  The  order,  title,  and  the 
notes  of  the  illustrations-  should  be  underneath  the  corresponding  frame 
and  translated  into  English.  The  drawing  paper  cannot  be  pasted  onto 
the  manuscript,  a  separate  pamphlet  of  hte  illustrations  shall  be  at¬ 
tached. 

4.  Equations  all  without  exception  shall  be  centered  on  the  man¬ 
uscript,  the  number  of  the  equation  shall  be  enclosed  in  (  )  on  the 
right  margin,  without  a  dotted  line. 

5.  Reference  literature  should  be  listed  by  number  in  back  of  th 
text,  according  to  the  exact  order  used  in  the  text;  only  list  primary 
and  published  literature.  It  shall  be  written  in  the  following  pattern 
(number)  Author,  Book  Title,  Publisher,  Year  Published,  Page  Numbers, 
(number).  Author,  Article  Subject  Title,  Periodical  Name,  Volume  Number 
Publication  Number,  Year,  Pages. 

6.  The  article  must  include  an  abstract  in  Chinese  and  English 
(number  of  words  not  to  exceed  300),  English  spelling  of  the  author's 
work  unit,  and  Chinese  spelling  in  pinyin  of  the  author's  name. 
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III.  Other 


1 .  If  submitted  manuscripts  have  been  read  at  a  certain  conference, 
have  appeared  in  another  journal,  or  have  already  been  published,  this 
must  be  made  clear. 

2.  After  submitted  manuscripts  have  been  examined,  all  essays 
that  are  carried  in  "Journal  of  Aeronautics"  will  receive  payment 
according  to  requlations. 

3.  Submitted  manuscripts  may  be  mailed  to:  Beijing,  Institute  Way, 
China  Aeronautics  Institute,  Journal  of  Aeronautics  Editorial  Depart¬ 
ment.  Be  sure  not  to  mail  it  to  an  individual. 
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